#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 #include #include 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 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: