MuonTimeModel.h

Go to the documentation of this file.

#ifndef _utl_MuonTimeModel_h_
#define _utl_MuonTimeModel_h_

#include <cmath>


namespace utl {
  class TabulatedFunction;
  class VRandomSampler;
  class RandomEngine;
}


namespace utl {

  class MuonTimeModel {

  public:
    MuonTimeModel(utl::RandomEngine& randomEngine,
                  const double theta = 0,
                  const double logE = 19,
                  const int primary = 1,
                  const bool flagAngularFactorDa_in = true);

    MuonTimeModel(utl::RandomEngine& randomEngine,
                  const double theta,
                  const utl::TabulatedFunction *LogzDist,
                  const bool flagAngularFactorDa_in = true,
                  const bool flagDecayFactor_in = true);

    ~MuonTimeModel();

    void Info();
    void SetCoordinates(const double r, const double zeta); // zeta rad
    void SetCoordinates(const double r, const double zeta, const double delta);

    double dNdz(const double z);
    double e_logt(const double logt);
    double g_logt(const double logt);
    double TotaldNdlogt(const double logt);
    double TotaldNdt(const double t);
    double GetDeltaTime();

    double GetFirstTime(const int n = 1);
    double GetTimes(const int n, double* const at);
    double GetMeanTime(const int n = 1);
    void GetFirstAndMeanTime(double& t_first, double& t_mean, const int n = 1);
    double GetLastTime(const int n = 1);
    void GetMeanAndRMSOfFirstTime(double& mean_t1,
                                  double& rms_t1,
                                  const int n = 1,
                                  const int stats = 1000);

  private:
    void UpdateModel();
    void DefaultSettings();

    double L_Z(const double z) const { return sqrt(z*z + fR*fR); }

    //This function was not in the original model
    //It was added to be able to use the PRODUCTION DISTRIBUTION h(z)
    //instead of the PRODUCTION DISTRIBUTION at ground dN/dz
    //depends on the flagDecayFactor
    double DecayFactor_Z(double Z);
    double cosaDa(double z);
    double Z_t(double t);
    double dZdt(double t);
    double dNdlogz(double logz);

    double g_t(double t);

    double E(double t, double x);
    double dEdt(double t, double x);
    double dNdE(double E, double x);
    double dNdlogE(double logE, double x);
    double e_t(double t);

    // *****************************************************************
    // L. Cazon parametric dN/dz model
    //
    double
    GetParametricLogMean(const double cosTheta, const double /*logE*/,
                         const int /*primary*/)
      const
    {
      return
        +5.682 + cosTheta*(
          -5.524 + cosTheta*(
            +9.180 + cosTheta*(
              -8.277 + cosTheta*(
                +2.797))));
    }

    double
    GetParametricLogSigma(const double cosTheta, const double /*logE*/,
                          const int /*primary*/)
      const
    {
      return
        +0.03266 + cosTheta*(
          +0.456 + cosTheta*(
            -1.084 + cosTheta*(
              +1.199 + cosTheta*(
                -0.4366))));
    }

    /*
      void MuonTimeModel::GetParametricNorm(const double cosTheta,
      const double logE,
      const double primary) {
      const double norm =
      -3.673E5
      +4.467E7 *cosTheta
      -9.47E7  *cosTheta*cosTheta
      +2.917E8 *cosTheta*cosTheta*cosTheta
      -1.781E8 *cosTheta*cosTheta*cosTheta*cosTheta;
      return norm;
      }
    */

    double
    GetParametricLogLambda(const double cosTheta, const double /*logE*/,
                           const int /*primary*/)
      const
    {
      const double logLambda =
        -0.01111 + cosTheta*(
          +0.2581 + cosTheta*(
            -0.7385 + cosTheta*(
              +3.362 + cosTheta*(
                -2.261))));
      if (logLambda < 0.03)
        return 0.03;
      return logLambda;
    }

    // *******************************************************************
    // Data members

    utl::RandomEngine& fRandomEngine;

    bool flagUserdNdlogz;
    int fNLogSteps;
    static const double fgM2;
    static const double fgPk;
    static const double fgKappa;
    static const double fgLambda;
    static const double fgQ;
    static const double fgGamma;

    // -- Basic shower parameters --
    double fTheta;
    double fCosTheta;
    double fLogE;
    int fPrimary;

    double fzMean;

    double fDelta;
    double fR;
    double fZeta;

    bool flagAngularFactorDa;
    bool flagDecayFactor;

    utl::VRandomSampler* fGRDGeometricalLogtDist;
    utl::VRandomSampler* fGRDKinematicalLogtDist;

    utl::TabulatedFunction* fGeometricalLogtDist;
    utl::TabulatedFunction* fKinematicalLogtDist;
    utl::TabulatedFunction* fLogzDist;

    double logstep_g;
    double logstep_e;
    double logt_g_up;
    double logt_g_low;
    double fLogt_e_up;
    double fLogt_e_low;

  private:
    // utility functions

    static void gser(double&, const double, const double, double&);
    static void gcf(double&, const double, const double, double&);
    static double GammaP(const double, const double);
    static double GammaQ(const double, const double);
    static double ErrFC(const double x);
    static double ErrF(const double);

  };

}


#endif