#ifndef _utl_PhysicalFuctions_h_
#define _utl_PhysicalFuctions_h_

/**
   \file
   Physical functions for common use in the framework
   
   \author Luis Prado Jr
   \version $Id: PhysicalFunctions.h 23499 2013-04-24 09:17:43Z roderigues $
   \date 19 May 2004
*/

#include <cmath>
#include <utl/Math.h>
#include <utl/PhysicalConstants.h>


namespace utl {
  
  /// Calculate the Gaisser-Hillas function
  double GaisserHillas(const double x, const double x0, 
                       const double xMax, const double nMax, 
                       const double lambda);

  /** Inverse of the Gaisser-Hillas function
      branch can be either -1 or +1
    \param n: the target particle number
    \param x0: slant depth of first interaction
    \param xMax: slant depth of shower maximum
    \param nMax: particle number at shower maximum
    \param lambda: radiation length
    \returns the X where GH(X) == n
   */
  template<int side>
  double
  GaisserHillasInverse(const double n, const double x0,
                       const double xMax, const double nMax,
                       const double lambda)
  {
    const double sxMax = (xMax - x0) / lambda;
    const double sn = n / nMax;
    const double x =
      -sxMax * LambertW<(side-1)/2>(-std::pow(sn, 1/sxMax) / M_E);
    return x * lambda + x0;
  }

  /// Calculate the electron energy for a relativistic beta
  inline double Energy(const double beta)
  { return kElectronMass / std::sqrt(1 - Sqr(beta)); }

  /// Calculate the electron energy versus the relativistic beta
  inline double Beta(const double energy)
  { return std::sqrt(1 - Sqr(kElectronMass/energy)); }

  /// Calculate the refraction index for a given depth
  double RefractionIndex(const double x);

  /// Wavelength-dependent index of refraction for humid air 
  double RefractionIndex(const double wl,
                         const double temperature,
                         const double pressure,
                         const double vaporPressure);

  /// Evaluate the saturation vapor pressure over ice or water
  double SaturationVaporPressure(const double temperature);
  
  /// Calculate the electron Cherenkov threshold energy for refraction index
  inline double CherenkovThreshold(const double nRef)
  { return kElectronMass / std::sqrt(1 - 1/Sqr(nRef)); }

  /// General definition of shower age
  inline double ShowerAge(const double slantDepth, const double showerMax)
  { return 3 / (1 + 2*showerMax/slantDepth); }

  /// The Moliere Radius (2 radiation length above obs-level, GAP-1998-002)
  double MoliereRadius(const double temperature,
                       const double pressure,
                       const double cosTheta = 0);

  /// parameter a of D. Gora et al., Astropart. Phys. 24 (2006), 484
  inline double GoraAParameter(const double age)
  { return (((5.151*age - 28.925)*age + 60.056)*age - 56.718)*age + 22.331; }

  /// parameter b of D. Gora et al., Astropart. Phys. 24 (2006), 484
  inline double GoraBParameter(const double age)
  { return (-1.039*age + 2.251)*age + 0.676; }

  /// CDF of lateral light distribution from Gora et al., APP24 (2006), 484
  /// input: distance to shower axis in units of moliere radii and shower age
  inline double GoraCDF(const double rStar, const double age) {
    const double a  = GoraAParameter(age);
    const double b  = GoraBParameter(age);
    return 1- pow(1+a*rStar,-b);
  }

  /// PDF of lateral light distribution from Gora et al., APP24 (2006), 484
  /// input: distance to shower axis in units of moliere radii and shower age
  inline double GoraPDF(const double rStar, const double age) {
    const double a  = GoraAParameter(age);
    const double b  = GoraBParameter(age);
    return b*a*pow(1+a*rStar,-b-1);
  }

  /// Fraction of electrons above energy cutoff enCut (in MeV) at age = 1
  double ElectronsAboveCut(const double enCut);
  
  /// Parametrization for the average energy deposit per particle
  double EnergyDeposit(const double age, const double enCut);
  
  /// Parametrization for the average energy deposit per particle
  inline double EnergyDeposit(const double depthX, const double xMax, const double enCut)
  { return EnergyDeposit(ShowerAge(depthX, xMax), enCut); }

  /// Model of the signal uncertainty (as given in GAP 2012-012)
  inline double SignalUncertaintyFactor(const double cosTheta)
  { return 0.34 + 0.46 / cosTheta; }

  enum EInvEInteractionModel {
    eQGSJET01, 
    eSIBYLL21,
    eGOLDEN
  };
    
  enum EInvECompositionModel {
    eMixedComposition, 
    eProton, 
    eIron,
    eGolden
  };

  /// invisible energy factor, finv=Etot/Eem, given Eem
  double
  InvisibleEnergyFactor(const double Eem, 
			const EInvEInteractionModel iMod, 
			const EInvECompositionModel iCompo);
 
  /// derivative of invisible energy factor dfinv/dEem given Eem
  double
  InvisibleEnergyFactorDerivative(const double Eem, 
				  const EInvEInteractionModel iMod, 
				  const EInvECompositionModel iCompo);
  
  /// variance of invisible energy factor due to shower-to-shower fluctuations
  double
  InvisibleEnergyFactorVariance(const double eTot);

}


#endif

// Configure (x)emacs for this file ...
// Local Variables:
// mode: c++
// compile-command: "make -C .. -k"
// End: