PhysicalFunctions.h

Go to the documentation of this file.

#ifndef _utl_PhysicalFuctions_h_
#define _utl_PhysicalFuctions_h_

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


namespace utl {

  double GaisserHillas(const double x, const double x0,
                       const double xMax, const double nMax,
                       const double lambda);


  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;
  }


  inline double Energy(const double beta)
  { return kElectronMass / std::sqrt(1 - Sqr(beta)); }


  inline double Beta(const double energy)
  { return std::sqrt(1 - Sqr(kElectronMass/energy)); }

  namespace RefractionIndex {

    double LorentzLorentz(const double verticalDepth);

    double GladstoneDale(const double density,
                         const double densityAtSeaLevel = kAirDensitySeaLevel,
                         const double refractiveIndexAtSeaLevel = kRefractiveIndexSeaLevel);

    double Ciddor95(const double wl,
                    const double temperature,
                    const double pressure,
                    const double vaporPressure);

  }


  double SaturationVaporPressure(const double temperature);


  inline double CherenkovThreshold(const double nRef)
  { return kElectronMass / std::sqrt(1 - 1/Sqr(nRef)); }


  inline double ShowerAge(const double slantDepth, const double showerMax)
  { return 3 / (1 + 2*showerMax/slantDepth); }


  double MoliereRadius(const double temperature,
                       const double pressure,
                       const double cosTheta = 0);


  double GoraCDF(const double rStar, const double age);


  double InverseGoraCDF(const double fraction, const double age);


  double GoraPDF(const double rStar, const double age);


  double ElectronsAboveCut(const double enCut);


  double EnergyDeposit(const double age, const double enCut);


  inline double EnergyDeposit(const double depthX, const double xMax, const double enCut)
  { return EnergyDeposit(ShowerAge(depthX, xMax), enCut); }




  bool
  QuadraticMaximumInterpolation(const std::vector<double>& x, const std::vector<double>& y,
                                double& xMax, double& yMax);

  namespace wcd {

    inline
    double
    PoissonFactor(const double sigmaFactor)
    {
      return std::max(1., 1/utl::Sqr(sigmaFactor));
    }


    enum ESignalVarianceModel {
      eGAP2012_012,
      eGAP2014_035
    };


    double
    SignalUncertaintyFactor(const ESignalVarianceModel model,
                            const double cosTheta);


    double
    SignalUncertainty(const ESignalVarianceModel model,
                      const double cosTheta, const double signal);


    double TriggerProbability(const bool totdMoPSEnabled,
                              const double lgExpectedSignal, const double sin2Theta);

  }


  namespace ssd {

    inline
    double
    PoissonFactor(const double sigmaFactor)
    {
      return 1/utl::Sqr(sigmaFactor);
    }

    enum ESignalVarianceModel {
      eGAP2018_025,
      eGAP2019_000 // until gap-note is published.
    };


    double
    SignalUncertaintyFactor(const ESignalVarianceModel model,
                            const double cosTheta);


    double
    SignalUncertainty(const ESignalVarianceModel model,
                      const double cosTheta, const double signal);

  }


  namespace InvisibleEnergy {

    enum EInteractionModel {
      eQGSJET01,
      eSIBYLL21,
      eDATA
    };


    enum ECompositionModel {
      eMixedComposition,
      eProton,
      eIron,
      eData
    };


    double
    EnergyLimits(const double energy);


    double
    ModelFactor(const EInteractionModel interaction,
                const ECompositionModel composition);


    const double*
    FitParameters(const ECompositionModel composition);


    double
    Factor(const double Eem,
           const EInteractionModel iMod,
           const ECompositionModel iCompo,
           const double CosTheta);


    double
    FactorDerivative(const double Eem,
                     const EInteractionModel iMod,
                     const ECompositionModel iCompo,
                     const double CosTheta);


    // ICRC 2013
    // variance of invisible energy factor due to shower-to-shower fluctuations
    // double FactorVariance(const double eTot);
    // ICRC 2019
    // variance of invisible energy factor due to shower-to-shower fluctuations and to the mass variance
    double
    FactorVariance(const double eCal, const double eTot);

  }


  // Enum specifying the high-energy hardronic interaction model to be used to estimate/predict Xmax.
  enum HadronicInteractionModel {
      eQGSJETII04,
      eEPOSLHC,
      eSibyll23d
  };


  namespace XmaxParam {
    /* Parameterization of the Xmax moments for CONEX simulated Xmax distributions.
       The equations from 10.1088/1475-7516/2013/02/026  
       The parameters (for the newest interaction models) from GAP2020_058. */
    
    double 
    Mean(const double energy, const double massNumber,
         const HadronicInteractionModel hadModel);

    double
    StandardDeviation(const double energy, const double massNumber, 
                      const HadronicInteractionModel hadModel);
  
  }
  
  
  namespace GumbelXmax {
    /* Parameters of the generalized gumbel function. 
       
       G(z) = 1 / sigma * lambda ^ lambda / Gamma(lambda) 
              * (exp(-lambda * z - lambda * exp(-z)))
       with z = (x - mu) / sigma
       
       The equations are from 10.1088/1475-7516/2013/07/050.
       The parameters (for the newest interaction models) 
       are from GAP2020_058. */

    // Does not correspond to the mean of the Xmax distribution
    double 
    Mu(const double energy, const double massNumber,
       const HadronicInteractionModel hadModel);

    // Does not correspond to the standard deviation of the Xmax distribution
    double
    Sigma(const double energy, const double massNumber, 
          const HadronicInteractionModel hadModel);

    double
    Lambda(const double energy, const double massNumber, 
           const HadronicInteractionModel hadModel);
  }

}


#endif