One way to do this goes something like this (for one methanol):
enter v1 vect xyz meoh 1 o meoh 1 h
traj firstu 101 nunit 1
mantime v1 normal
write v1 dumb time unit 200
* dummy title
The z-coordinate of v1 (4:th column written to unit 200) is cos(theta) for the OH of this methanol.
Note that if you used SHAKE during your simulation the O-H bond length is fixed to the value given in the parameter file (you can also measure it easily, by hand, using a molecular graphics program or the CHARMM commnand "Quick meoh 1 o meoh 1 h" (miscom.doc).