ClearAll[quietxAct];
SetAttributes[quietxAct,HoldAll];
xactoutfile=FileNameJoin[{NotebookDirectory[],"xact.out"}];
OpenWrite[xactoutfile];
Close[xactoutfile];
quietxAct[commands_]:=Block[{output=$Output},
$Output=OpenAppend[xactoutfile];
ReleaseHold@commands;
Close[$Output];
$Output=output;
]

quietxAct[
<<xAct`xTensor`;
<<xAct`xPert`;
<<xAct`xTras`;]

quietxAct[
DefConstantSymbol[dimd];
dimd=4;
DefManifold[spacetime,dimd,{μ,ν,ρ,σ,λ,α,β,γ,δ,τ,κ}];
DefMetric[-1,metric[-μ,-ν],CD,{";","𝒟"},SymCovDQ->True,PrintAs->"g"];
]

quietxAct[
DefMetricPerturbation[metric,hpert,ϵ];
PrintAs[hpert]^="(OverBar[h])";
Unprotect[IndexForm];
IndexForm[LI[X_]]:=ColorString[ToString[X],Red];
Protect[IndexForm];
LI/:delta[-LI[a_],LI[b_]]:=KroneckerDelta[a,b];
]

linsplit={hpert[LI[n_/;n>1],μ_,ν_]:>0};

quietxAct[

 DefTensor[h1[α,β],spacetime,Symmetric[{α,β}]];
 DefTensor[h2[α,β],spacetime,Symmetric[{α,β}]];
 
DefTensor[p1[μ],spacetime,PrintAs->"(Subscript[p, 1])"];
DefTensor[p2[μ],spacetime,PrintAs->"(Subscript[p, 2])"];

]

ΓEH=((Sqrt[-Detmetric[]]))(RicciScalarCD[])

$Post=ScreenDollarIndices

Γhh=Module[{res=ΓEH},
res=((Coefficient[Perturbed[res,2],ϵ,2])/(Sqrt[-Detmetric[]]));
res=res//ExpandPerturbation//ToFlat;
res=res/.linsplit//ToFlat//ContractMetric//ToCanonical//CollectTensors;
Return[res];
Return[res/.hpert->Zero];
];

(*In the following I am trying to do partial integration and discard total derivatives, to bring in terms like h*delh and delh*delh*)

resd=h1[μ,ν]*VarD[hpert[LI[1],μ,ν],CD][Γhh]//ContractMetric//ToCanonical//CollectTensors

resdd=h2[ρ,σ]*VarD[hpert[LI[1],ρ,σ],CD][resd]//ContractMetric//ToCanonical//CollectTensors