#ifndef _tls_MuonProfileUtilities_h_
#define _tls_MuonProfileUtilities_h_

#include <utl/AugerUnits.h>
#include <utl/Math.h>

namespace tls {

  /// kDeclination of geomagnetic field in Malargüe, 2007
  static const double kDeclination = 3.92*utl::degree;
  /// kInclination of geomagnetic field in Malargüe, 2007
  static const double kInclination = -35.09*utl::degree;

  /// Rotation angle needed to align the y-axis
  /// (perpendicular to the axis of the cascade) and transverse component
  /// of magnetic field. Original written by Ricardo Vazquez
  inline
  double
  GetBpsi(const double theta, const double phi)
  {
    using namespace std;
    //Magnetic field in our coordinate system.
    const double Bx =  sin(kDeclination)*cos(kInclination);
    const double By =  cos(kDeclination)*cos(kInclination);
    const double Bz = -sin(kInclination);

    //Rotate the magnetic field.
    const double Bx2 =  Bx*cos(phi) + By*sin(phi);
    const double By2 = -Bx*sin(phi) + By*cos(phi);
    const double Bz2 =  Bz;

    const double Bx3 = Bx2*cos(theta) - Bz2*sin(theta);
    const double By3 = By2;

    const double Bt = sqrt(Bx3*Bx3 + By3*By3);
    const double c_beta =  Bx3/Bt;
    const double s_beta =  By3/Bt;

    const double b_psi = atan2(s_beta,c_beta);

    return b_psi;
  }

  /// radial distance in shower front plane coordinate system
  inline
  double
  rPerpendicular(const double xpos, const double ypos,
                 const double theta, const double phi)
  {
    using namespace std;
    // align with shower azimuth
    double xRotate = xpos*cos(phi) + ypos*sin(phi);
    double yRotate = -xpos*sin(phi) + ypos*cos(phi);
    // project onto perpendicular plane
    double xPerpendicular = xRotate * cos(theta);
    double yPerpendicular = yRotate;

    double kRPerpendicular = sqrt(pow(xPerpendicular,2)+pow(yPerpendicular,2));

    return kRPerpendicular;
  }

  /// polar angle in shower front plane coordinate system
  inline
  double
  PsiPerpendicular(const double xpos, const double ypos,
                   const double theta, const double phi)
  {
    using namespace std;
    // align with shower azimuth
    double xRotate = xpos*cos(phi) + ypos*sin(phi);
    double yRotate = -xpos*sin(phi) + ypos*cos(phi);
    // project onto perpendicular plane
    double xPerpendicular = xRotate * cos(theta);
    double yPerpendicular = yRotate;

    return atan2(yPerpendicular,xPerpendicular);
  }

} // NS tls

#endif