首先先把每个连通块找出来。计算出每个连通块长度为 $t$ 的链的数量。

考虑简单问题：不对环的长度求和，只算方案数。

接下来考虑 $dp_{i,j}$ 表示看到第 $i$ 个环，目前整个的长度跟 $y$ 取 $\min$ 是 $j$ 的不同连接方案数。我们可以暂时不考虑顺序，在最后乘上 $\frac{cnt}{2}!$ 即可，其中 $cnt$ 为环数。

考虑转移，暴力转移就是对的。

证明：对于所有点数不超过 $\sqrt n$ 的连通块，其不同长度路径数至多为 $n\sqrt n$。而对于所有点数超过 $\sqrt n$ 的连通块，其不同长度路径数至多是 $y$。故总转移数是 $O((n+y)\sqrt n)$ 的。由于每次有 $y$ 个状态转移，总复杂度是 $O((n+y)y\sqrt n)$。

对于方案数，我们考虑 $f_{i,j}$，表示所有方案的环长度之和。每次 $f_{i,j}$ 既可以从 $f_{i-1,k}$ 转移，也可以从 $dp_{i-1,k}$ 通过长度的系数转移。复杂度一样。

#include <bits/stdc++.h>
#define int long long
#define double long double
#define lowbit(i) (i&(-i))
using namespace std;
const int mod=1e9+7,inv2=(mod+1)/2;
int qp(int a,int b){
	int ans=1;
	while(b){
		if(b&1) (ans*=a)%=mod;
		(a*=a)%=mod;
		b>>=1;
	}
	return ans;
}
int fac[1000005],inv[1000005];
void init(){
	fac[0]=1; for(int i=1;i<=1000000;i++) fac[i]=fac[i-1]*i%mod;
	inv[1000000]=qp(fac[1000000],mod-2); for(int i=999999;i>=0;i--) inv[i]=inv[i+1]*(i+1)%mod;
}
int C(int i,int j){
	if(i<0||j<0||i<j) return 0;
	return fac[i]*inv[j]%mod*inv[i-j]%mod;
}
vector<pair<int,int>> vc[100005];
int f[100005],siz[100005],cnt[3005][3005],val[3005][3005],rn[3005],cntt;
int dp2[3005][3005][2],y;
int find(int i){
	return f[i]==i?f[i]:f[i]=find(f[i]);
}
void dfs(int now,int fa,int rt,int dep){
	if(fa) cnt[rn[find(rt)]][min(dep,y)]++,(val[rn[find(rt)]][min(dep,y)]+=dep)%=mod;
	for(auto v:vc[now]){
		if(v.first==fa) continue;
		dfs(v.first,now,rt,dep+v.second);
	}
}
signed main(){
	init();
	int n,m,x,ans=0; cin>>n>>m>>x>>y;
	for(int i=1;i<=n;i++) f[i]=i;
	for(int i=1;i<=m;i++){
		int u,v,w; cin>>u>>v>>w;
		vc[u].push_back(make_pair(v,w));
		vc[v].push_back(make_pair(u,w));
		f[find(u)]=find(v);
	}
	for(int i=1;i<=n;i++) if(find(i)==i) rn[i]=++cntt;
	for(int i=1;i<=n;i++) dfs(i,0,i,0);
	{
		int st=min(y,(n-m)*x);
		dp2[0][st][0]=1;dp2[0][st][1]=(n-m)*x;
		for(int i=1;i<=cntt;i++){
			for(int k=0;k<=y;k++){
				if(cnt[i][k]){
					for(int j=0;j<=y;j++){
						(dp2[i][min(j+k,y)][0]+=dp2[i-1][j][0]*cnt[i][k])%=mod;
						(dp2[i][min(j+k,y)][1]+=dp2[i-1][j][1]*cnt[i][k]+dp2[i-1][j][0]*val[i][k])%=mod;
					}
				}
			}
		}
		(ans+=dp2[cntt][y][1]*fac[n-m-1]%mod*inv2)%=mod;
	}
	cout<<ans;
	return 0;
}