Thanks for explaining your motivation, because I think that the general problem as you stated it is impossibly hard, but that, fortunately, for the problem that you are really trying to tackle (the inverse problem in the calculus of variations), there is no need to solve this problem in this generality. If you are willing to take advantage of the intrinsic geometry of the spherized tangent bundle, you only need to solve a *much* easier problem. Indeed, if you don't take advantage of this geometry, solving the hard problem won't help for this purpose.

To see why the general problem is so hard, start with an orientable $2n$-manifold $M$, fix a volume form $\nu$ on $M$ (i.e., a nonvanishing $2n$-form on $M$), and let $\Omega = \pi_1^\ast(\nu)$ where $\pi_1:M\times S^1\to M$ is the projection on the first factor. If you did have manageable necessary and sufficient conditions to solve your problem as you have stated it, then, applying it to $\bigl(M\times S^1,\Omega\bigr)$, you'd have necessary and sufficient conditions for $M$ to be a symplectic manifold. Of course, this is known to be a very hard problem, and I don't know anyone who believes that we are going to have the solution to this general global problem anytime soon, even though, as you point out, the existence of local solutions is utterly trivial.

However, given that you are trying to solve the inverse problem in the calculus of variations, you don't really need to solve this general problem. You can take advantage of the fact that you are working on $\mathsf{R}M = \bigl(TM\setminus O_M)/\mathbb{R}^+$, which you call the 'spherized tangent bundle', but which I like to call the 'tangent ray bundle' of a smooth manifold $M^{n+1}$ (so I use '$\mathsf{R}$').

There are two geometric features of $\mathsf{R}M$ that are important: First, there is the foliation by the fibers ${\mathsf R}_x M$ for $x\in M$, each of which is diffeomorphic to the $n$-sphere, and the corresponding 'vertical' $n$-plane field $V\subset T(\mathsf{R}M)$, i.e., the kernel of the differential of the projection $\pi:\mathsf{R}M\to M$. Second, there is the $(n{+}1)$-plane field $C\subset T(\mathsf{R}M)$ that contains the vertical plane field $V$ and has the property that $\pi'(r)(C_r)= \mathbb{R}\cdot r \subset T_{\pi(r)}M$. The plane fields $V$ and $C$ are canonical in the sense that $V$ and $C$ are preserved under the natural action on $\mathsf{R}M$ of the diffeomorphisms of $M$ (and they are the only plane fields that are preserved). There is also an involution $\iota:\mathsf{R}M\to\mathsf{R}M$ that sends each ray to its opposite ray, but this won't be important in what follows.

Now, an *oriented path geometry* on $M$ is, by one definition, a choice of a line bundle $E\subset C$ over $\mathsf{R}M$ such that $C = E\oplus V$. The $2n$-parameter family of curves in $\mathsf{R}M$ that are tangent to $E$ (and hence foliate $\mathsf{R}M$) can be canonically oriented so that they then $\pi$-project to $M$ to be a $2n$-parameter family of oriented curves, one through each point oriented tangent to each ray based at that point. The so-called inverse problem in the calculus of variations is to determine whether there exists a (first-order) nondegenerate Lagrangian for oriented curves in $M$ such that the extremals of that Lagrangian are exactly the oriented curves generated by $E$. [A *(classical) path geometry* is simply an oriented path geometry that is invariant under $\iota$.]

In this context, a *first-order Lagrangian* is a section $\lambda$ of the line bundle $(C/V)^*\to \mathsf{R}M$. The reason is that, for any immersed curve $\gamma:[0,1]\to M$, its tangential lifting to $\mathsf{R}M$ defined by $\mathsf{R}\gamma(t)= \mathbb{R}^+\cdot\gamma'(t)$ then can be used to pull back $\lambda$ to $[0,1]$, so that it can be integrated over that interval, thus defining a first-order functional on oriented, immersed curves in $M$. (The point is that $(C_r/V_r)^\ast$ is naturally isomorphic to the dual of the tangent line $\mathbb{R}\cdot r$.)

How does one determine $E$ from a given $\lambda$? The process is as follows: First, one notes the classical lemma that, in this form, says that, for any given Lagrangian $\lambda$ on $\mathsf{R}M$, there exists a unique $1$-form $\delta\lambda$ on $\mathsf{R}M$ with the following properties: First, $V$ is in the kernel of $\delta\lambda$, so that it makes sense to evaluate $\delta\lambda$, which is a section of $(T/V)^\ast$, as an element of $(C/V)^\ast$; second, this evaluation agrees with $\lambda$; and, third, $\bigl(d(\delta\lambda)\bigr)(v,w)=0$ whenever $v,w\in C$. The $1$-form $\delta\lambda$ is known as the *Poincaré-Cartan form* of the Lagrangian $\lambda$. The mapping $\delta:C^\infty((C/V)^\ast)\to C^\infty((T/V)^\ast)$ is a linear, first-order operator.

We say that $\lambda$ is *nondegenerate* if $(d(\delta\lambda))^n$, which is a closed $2n$-form on $\mathsf{R}M$, is nonvanishing. In this case, since $\mathsf{R}M$ has dimension $2n{+}1$, there is a line bundle $E_\lambda\subset T(\mathsf{R}M)$ that is the kernel of $d(\delta\lambda)$, and it is not difficult to show that $E_\lambda\oplus V = C$ and that, moreover, the oriented curves defined by the oriented path geometry $E_\lambda$ are the extremal curves (with respect to compactly supported variations) of the functional defined by the Lagrangian $\lambda$.

So, the inverse problem can be understood as, given an oriented path geometry $E$, find those nondegenerate Lagrangians $\lambda$ such that $E= E_\lambda$.

Now, one can imagine a gradual attack on this problem that tries to swim back 'upstream' from $E$ to $\lambda$. First, note that if a nondegenerate $\lambda$ exists such that $E=E_\lambda$, then there will be a closed nonvanishing $2n$-form $\Omega = \bigl(d(\delta\lambda)\bigr)^n$ on $\mathsf{R}M$ whose kernel is $E$. Then one must write $\Omega$ in the form $\Omega = \omega^n$ for some closed $2$-form $\omega$. Of course, doing this in general is exactly the problem of the OP. Supposing this could be done, though, that's far from enough because there's no guarantee that $\omega$ would vanish when restricted to each $(n{+}1)$-plane $C_r$, and this would be necessary if $\omega = d(\delta\lambda)$. You'd have to go back and try a different $\omega$, or worse, consider all possible $\omega$s and try to select one that does do what you want. This seems hopelessly difficult.

However, it turns out that you don't really need to do any of this. There is a better way to proceed: First, write the identity map of $E$ in the form $ X\otimes \xi$ where the vector field $X$ on $\mathsf{R}M$ is a nonvanishing section of $E$ (this can always be done and $X$ will be unique up to a scalar multiple) and $\xi$ is the dual section of $E^\ast$. Define a linear, second-order differential operator
$$
D_E:C^\infty\bigl((C/V)^\ast\bigr)\to C^\infty\bigl((T/C)^\ast\otimes E^\ast\bigr)
$$
(the vector bundle $T/C$ has rank $n$ over $\mathsf{R}M$) by the rule
$$
D_E\lambda = i_X\bigl(d(\delta\lambda)\bigr) \otimes \xi\ .
$$
(It is not hard to verify that $D_E\lambda$ is well-defined and really is a section of $(T/C)^\ast\otimes E^\ast$, i.e., an $E^\ast$-valued $1$-form on $\mathsf{R}M$ that vanishes on elements of $C$.) Then, by construction, $E=E_\lambda$ if $\lambda$ is a nondegenerate Lagrangian that satisfies $D_E\lambda=0$.

The point is that the inverse problem is cast as a linear, second-order PDE for the unknown Lagrangian $\lambda$, with the extra condition that one is only interested in nondegenerate Lagrangians. Locally, it is $n$ equations for $1$ unknown.

This PDE system is determined when $n=1$ and is always locally solvable. However, when $n>1$, this is an overdetermined problem, and, for the generic $E$, there are no nondegenerate solutions to $D_E\lambda=0$. (I believe that it was Jesse Douglas in the 1930s who first did a serious, detailed study of this overdetermined problem, and he exhibited path geometries in the case $n=2$ for which there were no nondegenerate solutions to this equation.)

In general, the techniques of exterior differential systems provide methods for finding, or at least describing, the general solutions for a given specific $E$. Fortunately, these methods are much easier to apply and study than the problem posed by the OP.

There is, of course, an enormous literature on this subject (not always written in the best notation, I have to say). If there is interest, I can give a (very abbreviated) bibliography.