C++ Site::physT方法代码示例

/*!
  Non-Zero Momentum nucleon  from wall sources
*/
void Nuc2pt::calcWallMomNucleon(QPropW& quark, QPropWMomSrc& mom_quark)
{
  char *fname = "Nuc2pt()";
  ThreeMom sink_mom(mom_quark.Mom()) ;
  sink_mom.conj();
  //The momentum on the sink has to be complex conjugated!!
  // p --> - p 
  
  srcdesc[0]=srcdesc[1]=srcdesc[2]=srcdesc[3]=-1 ;
  srcdesc[4]=srcdesc[5]=-1 ; // TIZB
  switch (quark.SrcType())
    {
    case WALL:
      srcdesc[0] = quark.SourceTime() ;
      source="WALL" ;
      break ;
    default:
      ERR.General(cname,fname,"Bad source type\n") ;
    }
  
  quark_mass = quark.Mass() ;
    
  WilsonMatrix Q1;  
  WilsonMatrix Q2;
  WilsonMatrix Q3;  
  
  SnkType=POINT;
  sprintf(mom,"%i %i %i", mom_quark.Mom().cmp(0), 
	  mom_quark.Mom().cmp(1), mom_quark.Mom().cmp(2) );
  Site s ;
  for(s.Begin();s.End();s.nextSite())
    {
      int t = s.physT() ;
      Q1 =  quark[s.Index()] ;  
      if(HalfFerm) Q1.PParProjectSink() ;
      Q2 = Q1 ;
      // here is the sink so I need exp(ipx) 
      // I did conjugate the momentum at the construction
      Q3 = sink_mom.Fact(s)*mom_quark[s.Index()]  ; 
      if(HalfFerm) Q3.PParProjectSink() ;
      calc_nuc_(Q1, Q2, Q3,t) ;	  
    }
  plus_parity.GlobalSum() ;
  minus_parity.GlobalSum() ;
}

/*! 
  Computes
  \f[
   Trace(seqQ \, 
     \gamma_{\mu_1} \, \gamma_{\mu_2} \, \gamma_{\mu_3}\, 
     e^{-i mom . x}\, Quark)
  \f]
  where \f$\mu\f$ is the member int gamma and ThreeMom mom is
  the momentum injected at the operator. 
**/
void Nuc3ptGamma::InsertOp(CorrFunc& tmp,QPropW& seqQ, QPropW& Quark)
{
  char *fname = "InsertOp()";

  VRB.Func(cname, fname);

  Site s ;
  for(s.Begin();s.End();s.nextSite())
    {
      int t(s.physT()) ;
      WilsonMatrix sq(seqQ[s.Index()]) ;
      for(int mu(0);mu<G.N();mu++)
	sq.gr(G[mu]) ; // Multiply by gamma_mu 
      Complex cc(mom.Fact(s)) ; // The momentum factor exp(-ipx)
      cc *= Trace(sq,Quark[s.Index()]) ;
      tmp[t] += cc ;
    }
}

void AlgFourierProp::calcmom(const FourMom& p)
{
  char *fname="calcmom(const FourMom&)";
  VRB.Func(cname,fname);

  Float   p1( p.x() + 0.5*bc[0] );
  Float   p2( p.y() + 0.5*bc[1] );
  Float   p3( p.z() + 0.5*bc[2] );
  Float   p4( p.t() + 0.5*bc[3] );
 
  const Float PI(3.141592654);
 
  p1 *= 2.0*PI/(GJP.XnodeSites()*GJP.Xnodes());
  p2 *= 2.0*PI/(GJP.YnodeSites()*GJP.Ynodes());
  p3 *= 2.0*PI/(GJP.ZnodeSites()*GJP.Znodes());
  p4 *= 2.0*PI/(GJP.TnodeSites()*GJP.Tnodes());

  Site site;

  while ( site.LoopsOverNode() )
    {
      const Float px( site.physX()*p1 );
      const Float py( site.physY()*p2 );
      const Float pz( site.physZ()*p3 );
      const Float pt( site.physT()*p4 );
      
      const Float pdotx( px + py + pz + pt );


      Rcomplex fact( cos(pdotx), -sin(pdotx) );

      //--------------------------------------------
      //  Using AddMult is logically equivalent to
      //     momprop += fact * (*prop)[i]; ,
      //  but cuts down on the number of loops over
      //  the data. The time taken is reduced by
      //  a factor of two.
      //--------------------------------------------
      
      momprop.AddMult(fact,(*prop)[site.Index()]);
      
    }
}

/*!
  Calculate zero momentum nucleon two point functions
*/
void Nuc2pt::calcNucleon(QPropW& quark, int dohalf){
  char *cname = "Nuc2pt" ;
  char *fname = "calcNucleon(QPropW&)" ;

  strcpy(mom,"0 0 0");  

  HalfFerm=dohalf;

  srcdesc[0]=srcdesc[1]=srcdesc[2]=srcdesc[3]=-1 ;
  srcdesc[4]=srcdesc[5]=-1 ; // TIZB
  switch (quark.SrcType())
    {
    case WALL:
      srcdesc[0] = quark.SourceTime() ;
      source="WALL" ;
      break ;
    case BOX:
      srcdesc[0] = quark.BoxSrcStart() ; 
      srcdesc[1] = quark.BoxSrcEnd() ;
      srcdesc[2] = quark.SourceTime() ;
      if( quark. BoxSrcUseXYZOffset() ){ // use QPropWArg.{x,y,z} as offset
	srcdesc[3] = quark.PointSrcX() ; 
	srcdesc[4] = quark.PointSrcX() ; 
	srcdesc[5] = quark.PointSrcX() ; 
      }
      source="BOX" ;
      break ;
    case POINT:
      srcdesc[0] = quark.PointSrcX() ; 
      srcdesc[1] = quark.PointSrcY() ;
      srcdesc[2] = quark.PointSrcZ() ;
      srcdesc[3] = quark.SourceTime() ;
      source = "POINT" ;
      break ;
    case GAUSS_GAUGE_INV:
      srcdesc[0] = quark.Gauss_N() ; 
      srcdesc[1] = int(100*quark.Gauss_W()) ;
      srcdesc[2] = quark.SourceTime() ;
      source = "GAUSS" ;
      break ;
    default:
      ERR.General(cname,fname,"Bad source type\n") ;
    }
  quark_mass = quark.Mass() ;
    
  WilsonMatrix Q1;  
  WilsonMatrix Q2;
  WilsonMatrix Q3;
  WilsonMatrix Q4;
  WilsonMatrix Q5;
  
    switch (SnkType)
      {
      case POINT:
      case GAUSS_GAUGE_INV:	
	{
	  if(SnkType==GAUSS_GAUGE_INV)
	    sink="GAUSS" ;
	  else
	    sink="POINT" ;
	  Site s ;
	  for(s.Begin();s.End();s.nextSite())
	    {
	      int t = s.physT() ;
	      Q1 = quark[s.Index()] ;
	      if(HalfFerm) Q1.PParProjectSink() ;
	      Q2 = Q1 ; 
	      Q3 = Q1 ; 
	      Q4 = Q1 ; 
	      Q5 = Q1 ; 
	      calc_nuc_(Q1, Q2, Q3, Q4, Q5, t, dohalf) ;
	    }
	  plus_parity.GlobalSum() ;
	  minus_parity.GlobalSum() ;
	}
	break ;
      case WALL:
	{
	  sink="WALL" ;

	  if(!quark.GFixedSnk())
	    ERR.General(cname,fname,"Wall sink needs gauge fixing\n") ;

	  for(int t = 0 ; t<plus_parity.TimeSize();t++)
	    {
	      Q1 =  quark.WallSinkProp(t) ;
	      if(HalfFerm) Q1.PParProjectSink() ;
	      Q2 = Q1 ;
	      calc_nuc_(Q1, Q2, Q1, t) ;
	    }
	}
	break ;
      default:
	ERR.General(cname,fname,"Unknown Sink Type\n") ;
      } 
}

//.........这里部分代码省略.........
      aQ=seqQ[s.Index()] ;
      // Apply the gamma matrices on the sequential  propagator
      // We multiply the source indices of the sequential propagator
      // which are the anti-quark
      for(int mu(0);mu<G.N();mu++)
	aQ.gr(G[mu]) ; // Multiply by gamma_mu 

      // Loop over all derivative terms on the Quark
      for(D.Start();D.NotEnd();D.NextQuark()) 
	{
	  // Calculate the endpoints of the shift
	  D.CalcEndPoints(aq_vec,q_vec,s);
	  // The sink indices of the quark propagator are the quark
	  Q=Quark.GetMatrix(q_vec,buff) ;
	  
	 
	  D.CalcListOfLinks(dirs);

#ifdef DEBUG
printf("\n%s:: dirs ",cname);
for(int i(0);i<D.N();i++) printf("%i ",dirs[i]);
printf(" aq_vec: ");
for(int i(0);i<4;i++) printf("%i ",aq_vec[i]);
printf(" q_vec: ") ; 
for(int i(0);i<4;i++) printf("%i ",q_vec[i]);
#endif

	  // Start from the anti-quark and move along the path described
	  // by the list of dirs multiplying all the encountered links
	  // at the end of the day you should be at the position of the quark
          if(D.N()>1)  
	    {
	      U.ZeroMatrix();
	      lat.PathOrdProdPlus(U,aq_vec,dirs,D.N()) ;
	      //lat.PathOrdProd(U,aq_vec,dirs,D.N()) ;
	    }
	  else //if D.N()==1 PathOrdPlus does not work
	    {
	      int abs_dir = dirs[0]&3; //abs_dir is {0,1,2,3}
	      if(abs_dir == dirs[0]) // dir is forward
		U = *(lat.GetBufferedLink(aq_vec, abs_dir));
	      else // dir is backward. q_vec is now the site where the link is
		U.Dagger(*(lat.GetBufferedLink(q_vec, abs_dir)));
	    }

#ifdef DEBUG
Complex cboo(U.Tr()) ; 
printf("  TrU= (%g %g)\n",cboo.real(),cboo.imag());       
#endif

	  //Now that the gauge connection is constructed multiply it on
	  //the quark i.e the source indices of the anti-quark propagator
	  aQU.UMultSource(U,aQ) ;

	  Terms(D.DTermIndx(),s.Index()) = Trace(aQU,Q) ;

#ifdef DEBUG
printf("   (anti-Q, Q) aQU = (%6.3f %6.3f), (%6.3f %6.3f) (%6.3f %6.3f)\n",
       aQ.Trace().real(),aQ.Trace().imag(),
       Q.Trace().real(),Q.Trace().imag(),
       aQU.Trace().real(),aQU.Trace().imag());
printf("%s:: DTermIndx(): %i site: %i Value: ( %6.3f %6.3f ) \n",
       cname,D.DTermIndx(),s.Index(),
       Terms(D.DTermIndx(),s.Index()).real(), 
       Terms(D.DTermIndx(),s.Index()).imag());
#endif


	}
    }

  Complex tt ;
  // Loop over sites and construct the operator
  for(s.Begin();s.End();s.nextSite()) 
    {
      int t(s.physT()) ;
      for(D.Start();D.NotEnd();D.Next()) // Loop over all derivative terms
	{
	  //Calculate the endpoints of the shift
	  //aq_vec : anti quark location
	  // q_vec :      quark location
	  D.CalcEndPoints(aq_vec,q_vec,s);
	  Complex cc(Terms.GetTerm(D.DTermIndx(),aq_vec,tt)) ;

	  //Finaly the momentum factors
	  cc *= mom.Fact(s) ; // The momentum factor exp(-ipx)	
	  cc *= D.Fact()    ; /* The sign and the normalization factor 
				 for the current derivative term         */

	  tmp[t] += cc ;

#ifdef DEBUG
printf("%s:: DTermIndx(): %i site: %i Fact: %9.6f Value: ( %6.3f %6.3f ) \n",
       cname, D.DTermIndx(),s.Index(),D.Fact(),cc.real(),cc.imag());
#endif
	
	}
    }
  delete [] dirs ;
}