Based on the answer of Robert Bryant, I was able to derive some more conditions which look a bit more geometric. As above $\omega$ is the contact form, $\Delta$ the contact structure, $X$ is the reeb vector field and $\theta$ is the other 1-form which defines the plane distribution $P$. I will assume to start with that $P$ and $\Delta$ are distinct everywhere. I will work local coordinates. So they intersect in a one dimensional line bundle spanned by some vector field $Y$ so that $\theta(Y) = \eta(Y)=0$.

I derived two conditions. 

First one:

Now I denote by $\eta_f = f\eta$ and its Reeb vector field $R_f$. Now since $d\eta_f(R_f,\cdot)=0$, one has that $\theta(R_f)=0$ if and only if 
$d\eta_f = \theta \wedge \beta_f$ for some everywhere non-vanishing 1-form $\beta_f$. This means $d\eta_f|_P=0$. Taking $Y$ as above and $X$ any other vectorfield, expanding $d\eta_f = df \wedge \eta + f d\eta$ and evaluating this 2-form on $(X,Y)$ one gets the equation:

$$\mathcal{L}_Y(ln(f)) = \eta(\mathcal{L}_Y(\frac{X}{\eta(X)})) $$

This is the first equation, I don't know if it admits a global solution.

Second one: 

We know that $\theta \wedge d\eta_f=0$ by above. We let $Z$ be a vectorfield inside $\Delta$ distinct from $Y$. Then we evaluate this on $(Z,Y,R)$ to get

$$\mathcal{L}_Y(ln(f)) = \frac{\theta(R)}{\theta(X)}d\eta(Z,Y)$$

In this case the quantity $d\eta(Z,Y)$ can be related to how the contact structure $\Delta$ twists along the integral curves of $Y$ and I have hopes somehow that right hand side can be written as the differential of a function but I am not sure because here the ugly term is $\theta(R)$ which may vanish at some points.