MuonTimeModel.cc

Go to the documentation of this file.

#include "MuonTimeModel.h"

#include <utl/RandomEngine.h>
#include <utl/TabulatedFunction.h>
#include <utl/RandomSamplerFromPDF.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
#include <utl/Math.h>

#include <iostream>
#include <iomanip>

using namespace utl;
using namespace std;


/* Model Constants */
/* todo
   - add many comments (especially to constants)
   - only use meaningfull names or variables/functins
*/
const double MuonTimeModel::fgM2 = 0.011;
const double MuonTimeModel::fgPk = 0.0002;
const double MuonTimeModel::fgKappa = 0.8;
const double MuonTimeModel::fgLambda = 1;
const double MuonTimeModel::fgQ = 0.17;
const double MuonTimeModel::fgGamma = 2.6;


MuonTimeModel::MuonTimeModel(utl::RandomEngine& randomEngine,
                             const double theta,
                             const double logE,
                             const int primary,
                             const bool flagAngularFactorDa_in) :
  fRandomEngine(randomEngine),
  fTheta(theta),
  fCosTheta(cos(theta)),
  fLogE(logE),
  fPrimary(primary),
  flagAngularFactorDa(flagAngularFactorDa_in),
  flagDecayFactor(false),
  fGRDGeometricalLogtDist(0),
  fGRDKinematicalLogtDist(0),
  fGeometricalLogtDist(0),
  fKinematicalLogtDist(0),
  fLogzDist(0)
{
  flagUserdNdlogz = false;

  fzMean = pow(10, GetParametricLogMean(fCosTheta, fLogE, fPrimary));
  DefaultSettings();
}


MuonTimeModel::MuonTimeModel(utl::RandomEngine& randomEngine,
                             const double theta,
                             const TabulatedFunction* LogzDist,
                             const bool flagAngularFactorDa_in,
                             const bool flagDecayFactor_in) :
  fRandomEngine(randomEngine),
  fTheta(theta),
  fCosTheta(cos(theta)),
  flagAngularFactorDa(flagAngularFactorDa_in),
  flagDecayFactor(flagDecayFactor_in),
  fGRDGeometricalLogtDist(0),
  fGRDKinematicalLogtDist(0),
  fGeometricalLogtDist(0),
  fKinematicalLogtDist(0),
  fLogzDist(new TabulatedFunction(*LogzDist))
{
  flagUserdNdlogz = true;

  double sumXW = 0;
  double sumW = 0;
  for (unsigned int iBin = 0; iBin < fLogzDist->GetNPoints(); ++iBin) {
    sumXW += (*fLogzDist)[iBin].Y() * (*fLogzDist)[iBin].X();
    sumW += (*fLogzDist)[iBin].Y();
  }
  fzMean = pow(10, sumXW/sumW);

  DefaultSettings();
}


void
MuonTimeModel::DefaultSettings()
{
  fNLogSteps = 1000;
  fLogt_e_up = 5; // 3.23;
  fLogt_e_low = -6;
  fR = 1000;
  fZeta = 0;
}


MuonTimeModel::~MuonTimeModel()
{
  delete fGRDGeometricalLogtDist;
  delete fGRDKinematicalLogtDist;
  delete fGeometricalLogtDist;
  delete fKinematicalLogtDist;
  delete fLogzDist;
}


void
MuonTimeModel::Info()
{
  cout.setf(ios_base::fixed);
  cout << "\n"
          "*******************************\n"
          "* MODEL CONSTANTS             *\n"
          "*******************************\n"
          "* pk             "
       << setw(5) << setprecision(2)
       << fgPk*1000 << " GeV/km *\n"
          "* kappa         "
       << setw(5) << setprecision(1)
       << fgKappa << "         *\n"
          "* lambda        "
       << setw(5) << setprecision(1)
       << fgLambda << "         *\n"
          "* gamma         "
       << setw(5) << setprecision(1)
       << fgGamma << "         *\n"
          "* ptc or Q       "
       << setw(5) << setprecision(2)
       << fgQ << " GeV    *\n"
          "*******************************\n"
          "\n"
          "*******************************\n"
          "* SHOWER CONSTANTS            *\n"
          "*******************************\n"
          "* fTheta        "
       << setw(5) << setprecision(1)
       << fTheta/deg << "  deg    *\n"
          "*******************************\n"
          "\n"
          "*******************************\n"
          "* dN/d log10(z/m)             *\n"
          "*******************************\n";
  if (!flagUserdNdlogz) {
    const double cosTheta = cos(fTheta);
    cout << "* logMean /m      "
         << setw(5) << setprecision(3)
         << GetParametricLogMean(cosTheta, fLogE, fPrimary)
         << "       *\n"
            "* logSigma /m     "
         << setw(5) << setprecision(3)
         << GetParametricLogSigma(cosTheta, fLogE, fPrimary)
         << "       *\n"
            "* logLambda /m    "
         << setw(5) << setprecision(3)
         << GetParametricLogLambda(cosTheta, fLogE, fPrimary)
         << "       *\n";
  } else {
    cout << "* User Defined function       *\n"
            "* logMean /m      "
         << setw(5) << setprecision(3) << log10(fzMean) << "       *\n";
  }
  cout << "*******************************\n";
  if (flagAngularFactorDa)
    cout << "* D(a)=1                      *\n";
  else
    cout << "* D(a)=cos(a)+Delta/L         *\n";
  cout << "*******************************\n"
          "\n"
          "*******************************\n"
          "* COORDINATES                 *\n"
          "*******************************\n"
          "*fR          "
       << setw(10) << setprecision(1)
       << fR << "  m    *\n"
          "* fZeta         "
       << setw(5) << setprecision(1)
       << fZeta/deg << "  deg    *\n"
          "*******************************" << endl;
}


//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
MuonTimeModel::DecayFactor_Z(const double Z)
{
  if (flagDecayFactor)
    return pow(Z*Z + fR*fR, 0.5*(1 - fgGamma));
  else
    return 1;
}


double
MuonTimeModel::cosaDa(const double Z)
{
  if (Z < 0)
    return 0; //  muons back!?
  const double L = L_Z(Z);
  const double cosa2 = 1 - (fR*fR) / (L*L);

  if (flagAngularFactorDa)
    return sqrt(cosa2);
  return cosa2 + fDelta/L*sqrt(cosa2);
}


double
MuonTimeModel::Z_t(const double t)
{
  return 0.5*(fR*fR/(kSpeedOfLight*t) - kSpeedOfLight*t);
}


double
MuonTimeModel::dZdt(const double t)
{
  return 0.5*(fR*fR/(kSpeedOfLight*t*t) + kSpeedOfLight);
}


double
MuonTimeModel::dNdlogz(const double logz)
{
  if (flagUserdNdlogz) {
    if (logz < fLogzDist->XFront() ||
        logz > fLogzDist->XBack()-1)
      return 0;
    return fLogzDist->Y(logz);
  }

  const double logSigma = GetParametricLogSigma(fCosTheta, fLogE, fPrimary);
  const double logLambda = GetParametricLogLambda(fCosTheta, fLogE, fPrimary);
  const double logMean = GetParametricLogMean(fCosTheta, fLogE, fPrimary);

  if (logLambda <= 0) {
    const double e = (logz - logMean) / logSigma;
    static double sqrt2pi = 2.50662827463100024;
    return exp(-0.5*e*e) / (logSigma*sqrt2pi);
  }

  const double norm = 1;
  const double s2 = sqrt(2.);

  const double f1 = (logz - logMean) / logLambda;
  const double f2 = 0.5 * (logSigma*logSigma) / (logLambda*logLambda);
  const double f3 = (logz - logMean) / (s2*logSigma);
  const double f4 = logSigma / (s2*logLambda);

  return norm * exp(f1 + f2) * ErrFC(f3 + f4) / (2*logLambda);
}


double
MuonTimeModel::dNdz(const double z)
{
  if (z < 0)
    return 0;
  return 1 / (z*log(10.)) * dNdlogz(log10(z));
}


double
MuonTimeModel::g_t(const double t)
{
  const double Z = Z_t(t);
  if (t <= 0)
    return 0;
  return DecayFactor_Z(Z) * cosaDa(Z) * dNdz(Z+fDelta) * dZdt(t);
  //dZdt = dzdt
}


double
MuonTimeModel::g_logt(const double logt)
{
  const double t = pow(10, logt);
  return t * log(10.) * g_t(t);
}



double
MuonTimeModel::E(const double t, const double x)
{
  return 0.5 * fgPk * x * (1 + sqrt(1 + 2*fgM2/(fgPk*fgPk*x*kSpeedOfLight*t)));
}


double
MuonTimeModel::dEdt(const double t, const double x)
{
  const double pk2 = fgPk * fgPk;
  return 0.5 * fgPk * x /
    sqrt(1 + 2*fgM2 / (pk2 * x * kSpeedOfLight * t)) *
    fgM2 /
    (pk2 * x * kSpeedOfLight * t*t);
}


double
MuonTimeModel::dNdE(const double E, const double x)
{
  if (E < fgPk*x || x <= 0 || fR >= x)
    return 0;

  return (1 - fR*fR/(x*x)) * pow(fR/x, fgLambda) *
    pow(E, -fgGamma - fgKappa + fgLambda + 1) *
    pow(E - fgPk*x, fgKappa) *
    exp(-E*fR/(x*fgQ));
}


double
MuonTimeModel::dNdlogE(const double logE, const double x)
{
  const double E = pow(10, logE);
  return E * log(10.) * dNdE(E, x);
}


double
MuonTimeModel::e_t(const double t)
{
  if (t <= 0)
    return 0;

  const double x = L_Z(fzMean - fDelta);
  return dNdE(E(t, x), x) * dEdt(t, x);
}


double
MuonTimeModel::e_logt(const double logt)
{
  const double t = pow(10, logt);
  return t * log(10.) * e_t(t);
}


void
MuonTimeModel::SetCoordinates(const double r,
                              const double zeta)
{
  fR = r;
  fZeta = zeta;
  fDelta = fR * cos(fZeta) * tan(fTheta);
  UpdateModel();
}


void
MuonTimeModel::SetCoordinates(const double r,
                              const double zeta,
                              const double delta)
{
  fR = r;
  fZeta = zeta;
  fDelta = delta;
  UpdateModel();
}


void
MuonTimeModel::UpdateModel()
{
  const double LZ = (fzMean - fDelta);

  logt_g_up = log10(8000*(sqrt(fR*fR + LZ*LZ) - LZ));
  logt_g_low = log10(0.2*(sqrt(fR*fR + LZ*LZ) - LZ));

  logstep_g = (logt_g_up - logt_g_low) / fNLogSteps;
  logstep_e = (fLogt_e_up - fLogt_e_low) / fNLogSteps;

  delete fGeometricalLogtDist;
  fGeometricalLogtDist = new TabulatedFunction();
  delete fKinematicalLogtDist;
  fKinematicalLogtDist = new TabulatedFunction();

  for (int i = 0; i < fNLogSteps; ++i) {

    const double logtg = logt_g_low + logstep_g*i;
    const double logte = fLogt_e_low + logstep_e*i;

    fGeometricalLogtDist->PushBack(logtg, g_logt(logtg));
    fKinematicalLogtDist->PushBack(logte, e_logt(logte));

  }

  //------------------------

  delete fGRDGeometricalLogtDist;
  fGRDGeometricalLogtDist =
    new RandomSamplerFromPDF(*fGeometricalLogtDist, RandomSamplerFromPDF::eLinear);

  delete fGRDKinematicalLogtDist;
  fGRDKinematicalLogtDist =
    new RandomSamplerFromPDF(*fKinematicalLogtDist, RandomSamplerFromPDF::eLinear);
}


double
MuonTimeModel::GetFirstTime(const int n/*=1*/)
{
  double t_first = 1e25;
  for (int i = 0; i < n; ++i) {

    const double t_i =
      pow(10, fGRDGeometricalLogtDist->shoot(fRandomEngine.GetEngine())) +
      pow(10, fGRDKinematicalLogtDist->shoot(fRandomEngine.GetEngine()));

    if (t_i < t_first)
      t_first = t_i;
  }

  return t_first - fDelta/kSpeedOfLight;
}


double
MuonTimeModel::GetTimes(const int n, double* const at)
{
  double t_first = 1e25;
  for (int i = 0; i < n; ++i) {

    const double t_i =
      pow(10, fGRDGeometricalLogtDist->shoot(fRandomEngine.GetEngine())) +
      pow(10, fGRDKinematicalLogtDist->shoot(fRandomEngine.GetEngine()));

    at[i] = t_i - fDelta/kSpeedOfLight;
    if (t_i < t_first)
      t_first = t_i;
  }

  return t_first - fDelta/kSpeedOfLight;
}


double
MuonTimeModel::GetMeanTime(const int n/*=1*/)
{
  double t_mean = 0;
  for (int i = 0; i < n; ++i) {

    t_mean +=
      pow(10, fGRDGeometricalLogtDist->shoot(fRandomEngine.GetEngine())) +
      pow(10, fGRDKinematicalLogtDist->shoot(fRandomEngine.GetEngine()));

  }
  return t_mean/n - fDelta/kSpeedOfLight;
}


void
MuonTimeModel::GetFirstAndMeanTime(double& t_first,
                                   double& t_mean,
                                   const int n/*=1*/)
{
  t_mean = 0;
  t_first = 1e25;
  for (int i = 0; i < n; ++i) {

    const double t_i =
      pow(10, fGRDGeometricalLogtDist->shoot(fRandomEngine.GetEngine())) +
      pow(10, fGRDKinematicalLogtDist->shoot(fRandomEngine.GetEngine()));

    if (t_i < t_first)
      t_first = t_i;

    t_mean += t_i;
  }

  t_first -= fDelta / kSpeedOfLight;
  t_mean = t_mean/n - fDelta/kSpeedOfLight;
}


double
MuonTimeModel::GetLastTime(const int n/*=1*/)
{
  double t_last = 0;
  for (int i = 0; i <n; ++i) {

    const double t_i =
      pow(10, fGRDGeometricalLogtDist->shoot(fRandomEngine.GetEngine())) +
      pow(10, fGRDKinematicalLogtDist->shoot(fRandomEngine.GetEngine()));

    if (t_i > t_last)
      t_last = t_i;

  }
  return t_last - fDelta/kSpeedOfLight;
}


void
MuonTimeModel::GetMeanAndRMSOfFirstTime(double& mean_t1,
                                        double& rms_t1,
                                        const int n/*=1*/,
                                        const int stats/*=1000*/)
{
  double t1sum = 0;
  double t1sum2 = 0;

  for (int i = 0; i < stats; ++i) {
    const double t1 = GetFirstTime(n);
    t1sum += t1;
    t1sum2 += t1*t1;
  }

  mean_t1 = t1sum / stats;
  rms_t1 = sqrt(t1sum2/stats - mean_t1*mean_t1);
}


double
MuonTimeModel::GetDeltaTime()
{
  return fDelta / kSpeedOfLight;
}



double
MuonTimeModel::TotaldNdt(const double t)
{
  const int n = 850;
  const double tmin = 0;
  double tmax = 1700;

  if (t < tmax)
    tmax = t;

  const double h = (tmax - tmin) / n;
  double sum = 0;
  double t_i = tmin + 0.5*h;

  for (int i = 0; i < n; ++i) {
    sum += g_t(t - t_i) * e_t(t_i);
    t_i += h;
  }

  return sum * h;
}


double
MuonTimeModel::TotaldNdlogt(const double logt)
{
  const double t = pow(10, logt);
  return t * log(10.) * TotaldNdt(t);
}


/**************************************************
 *
 * Returns the incomplete gamma function P (a, x)
 * evaluated by its series representation as gamser.
 * Also returns ln G(a) as gln.
 *
 *  taken from Numerical Recipes in C
 *
 ***************************************************/

const int kITMAX = 100;      // max iterations
const double kEPS = 3e-7;    // relative accuracy
const double kFPMIN = 1e-30; // Number near the smallest representable floating-point number.


void
MuonTimeModel::gser(double& gamser, const double a, const double x, double& gln)
{
  gln = LogGamma(a);

  if (x <= 0) {
    if (x < 0)
      cerr << "x less than 0 in routine gser" << endl;
    gamser = 0;
  } else {
    double ap = a;
    double del = 1 / a;
    double sum = del;
    for (int n = 1; n <= kITMAX; ++n) {
      ++ap;
      del *= x/ap;
      sum += del;
      if (fabs(del) < fabs(sum)*kEPS) {
        gamser = sum * exp(-x + a*log(x) - gln);
        return;
      }
    }
    cerr << "a too large, ITMAX too small in routine gser" << endl;
  }
}


/**************************************************
 *
 * Returns the incomplete gamma function Q(a, x)
 * evaluated by its continued fraction represen-
 * tation as gammcf. Also returns ln (a) as gln.
 *
 *  taken from Numerical Recipes in C
 *
 ***************************************************/
void
MuonTimeModel::gcf(double& gammcf, const double a, const double x, double& gln)
{
  gln = LogGamma(a);
  //gammln(a);         // Set up for evaluating continued fraction

  double b = x + 1 - a;  // by modified Lentz's method (§5.2)
  double c = 1. / kFPMIN; // with b0 = 0.
  double d = 1 / b;
  double h = d;          //  Iterate to convergence.

  int i;
  for (i = 1; i <= kITMAX; ++i) {
    const double an = i * (a - i);
    b += 2;
    d = an*d + b;
    if (fabs(d) < kFPMIN)
      d = kFPMIN;
    c = b + an/c;
    if (fabs(c) < kFPMIN)
      c = kFPMIN;
    d = 1 / d;
    const double del = d*c;
    h *= del;
    if (fabs(del - 1) < kEPS)
      break;
  }
  if (i > kITMAX)
    cerr << "a too large, ITMAX too small in gcf" << endl;

  gammcf = exp(-x + a*log(x) - gln) * h; // Put factors in front.
}


/**************************************************
 *
 * Returns the incomplete gamma function P (a, x).
 *
 *  taken from Numerical Recipes in C
 *
 ***************************************************/
double
MuonTimeModel::GammaP(double a, double x)
{
  if (x < 0 || a <= 0)
    cerr << "gammp: Invalid arguments in x,a="
         << x << ' ' << a << endl;

  double g;
  double gln;
  if (x < a+1) { //   Use the series representation.
    gser(g, a, x, gln);
    return g;
  } else { //  Use the continued fraction representation
    gcf(g, a, x, gln);
    return 1 - g; //    and take its complement.
  }
}


/**************************************************
 *
 * Returns the incomplete gamma function Q(a, x)  1 - P (a, x).
 *
 *  taken from Numerical Recipes in C
 *
 ***************************************************/
double
MuonTimeModel::GammaQ(const double a, const double x)
{
  if (x < 0 || a <= 0)
    cerr << "Invalid arguments in routine gammq" << endl;

  double g;
  double gln;
  if (x < a+1) { // Use the series representation
    gser(g, a, x, gln);
    return 1 - g; // and take its complement.
  } else { //  Use the continued fraction representation.
    gcf(g, a, x, gln);
    return g;
  }
}


double
MuonTimeModel::ErrF(const double x)
{
  return x < 0 ? -GammaP(0.5, x*x) : GammaP(0.5, x*x);
}


double
MuonTimeModel::ErrFC(const double x)
{
  return x < 0 ? 1 + GammaP(0.5, x*x) : GammaQ(0.5, x*x);
}