#include <math.h>

double target_spin_angle( double z, double nominal ) {

   // Routine to calculate the spin angle as a function of the Z coordinate
   // based on the measured profile and a given "nominal" holding angle of
   // approximately 32 or 47 degrees
   //
   // NOTE: Analyses which strive to determine or verify the "nominal" angle
   // (like T20, e'p) should fit the data by varying the value for "nominal".
   //
   // This code allows a range of +/- 7 degrees about nominal which is
   // hopefully sufficient for searches.  Beyond that the nominal value is 
   // simply returned.
   //
   // Once an optimal value is found all other analyses should use that.
   //
   // Usage : z = event vertex in [cm]
   //         nominal = value close to 32.0 or 45.0 [deg]
   //         target_spin_angle is returned in degrees.
   //
   // Profiles are based on the 2005 and 2006 magnetic field maps done by Genya
   // using the Hall probe and later verified using the short compass.
   //
   // The fits were restricted to z values of +/- 20 cm.  Beyond this the
   // measurements go rapidly to very large angles.  Since most data is from
   // the central region of the target and cuts of +/- 20 cm are typical this
   // is probably okay but users should realise that if z input is beyond this
   // range the code just returns the nominal angle as the fits are no longer
   // valid and the behavior is undefined beyond this range.
   //
   //                                             - D.K. Hasell 19/4/06
   //

   double zz = z /20.0;  // Scale the z value so coefficients are reasonable.

   // If z is outside valid range for fits just return the nominal value.

   if( fabs( zz ) > 1.0 ) return nominal;

   // Choose the appropriate profile.

   if( fabs( nominal - 47.0 ) < 7.0 ) {

      // For a nominal angle near 47 degrees - 9th order polynomial.

      return nominal + zz *
	 ( -2.512581 + zz *
	   ( -29.581168 + zz *
	     ( 6.150861 + zz *
	       ( 93.978928 + zz *
		 ( -15.194133 + zz *
		   ( -104.159140 + zz *
		     ( 15.573366 + zz *
		       ( 50.518515 + zz * -5.611146 ) ) ) ) ) ) ) );
   }

   else if( fabs( nominal - 32.0 ) < 7.0 ) {

      // For a nominal angle near 32 degrees - 9th order polynomial.

      return nominal + zz *
	 ( -2.439392 + zz *
	   ( -44.961079 + zz *
	     ( 8.645792 + zz *
	       ( 134.282715 + zz *
		 ( -32.205690 + zz *
		   ( -160.727705 + zz *
		     ( 44.282327 + zz *
		       ( 77.849926 + zz *  -21.214353 ) ) ) ) ) ) ) );
   }

   else {

      // Nominal value out of range, simply return it.

      return nominal;

   }

   // That's All Folks !

}