UnivTimeKG.cc

Go to the documentation of this file.

#include "UnivTimeKG.h"

#include <utl/Math.h>
#include <utl/AugerUnits.h>

#include <Math/PdfFuncMathCore.h>
#include <Math/ProbFuncMathCore.h>
#include <Math/QuantFuncMathCore.h>

#include <fstream>
#include <cmath>


/*
  Implementation of the time model (base class)

  ===========
  parameters:
  ===========

  icomp: signal component number:
    pure muonic (0)
    pure electromagnetic (1)
    electromagnetic from muon decay (2)
    electromagnetic from hadron jets (3)

  r: distance to shower axis in shower plane CS
  psi: azimuth to shower axis in shower plane CS
  DX: distance to Xmax, parallel to shower axis
  t: relative time for certain bin in trace referred to plane front arriving at station

  Values are given in Offline units
*/


using namespace std;


namespace UnivTimeKG {


  unsigned int
  FindIndex(const vector<double>& array, const double value)
  {
//#warning this should be replaced with std::upper_bound or std::lower_bound
    unsigned int imin = 0;
    unsigned int imax = array.size() - 1;

    while (imax - imin > 1) {
      const unsigned int i = (imax + imin) / 2;
      if (array[i] <= value)
        imin = i;
      else
        imax = i;
      if (array[i] == value)
        break;
    }
    return imin;
  }


  double
  int1_2p(const double x0, const double y0, const double x1, const double y1, const double x)
  {
    return y0 + (y1 - y0) * (x - x0) / (x1 - x0);
  }


  double
  int1(const vector<double>& xs, const vector<double>& ys, const double x)
  {
    const unsigned int i = FindIndex(xs, x);
    const double x0 = xs[i];
    const double x1 = xs[i + 1];
    const double y0 = ys[i];
    const double y1 = ys[i + 1];
    return int1_2p(x0, y0, x1, y1, x);
  }


  double
  int1(const double xs[], const double ys[], const double x, const int n)
  {
    const vector<double> vxs(xs, xs + n);
    const vector<double> vys(ys, ys + n);
    return int1(vxs, vys, x);
  }


  double
  int2(const vector<double>& x1s, const vector<double>& x2s, const vector<vector<double> >& Ys, const double x1, const double x2)
  {
    const unsigned int j = FindIndex(x2s, x2);
    const double x21 = x2s[j];
    const double x22 = x2s[j + 1];

    const double fx1 = int1(x1s, Ys[j],   x1);
    const double fx2 = int1(x1s, Ys[j + 1], x1);
    const double fxy = int1_2p(x21, fx1, x22, fx2, x2);

    return fxy;
  }


  TimeModel::TimeModel()
  { }


  TimeModel::~TimeModel()
  { }


  void
  TimeModel::addInterpolationTable(const vector<double>& rs, const vector<double>& psis, const vector<vector<double> >& ys)
  {
    IntTable t;
    t.rs = rs;
    t.psis = psis;
    t.ys = ys;
    tables.push_back(t);
  }


  void
  TimeModel::addInterpolationTable(const std::string& filename)
  {
    ifstream file(filename.c_str());

    vector<double> rs;
    vector<double> psis;
    vector<vector<double> > ys;

    unsigned int Nr, Npsi;
    if (!(file >> Nr >> Npsi)) {
      cerr << "error reading file!" << endl;
      return;
    }

    double y = 0;

    for (unsigned int i = 0; i < Nr; ++i) {
      if (file >> y)
        rs.push_back(y);
      else {
        cerr << "error reading file!" << endl;
        return;
      }
    }

    for (unsigned int j = 0; j < Npsi; ++j) {
      if (file >> y)
        psis.push_back(y);
      else {
        cerr << "error reading file!" << endl;
        return;
      }
    }

    for (unsigned int i = 0; i < Nr; ++i) {
      vector<double> line;
      for (unsigned int j = 0; j < Npsi; ++j) {
        if (file >> y)
          line.push_back(y);
        else {
          cerr << "error reading file!" << endl;
          return;
        }
      }
      ys.push_back(line);
    }

    addInterpolationTable(rs, psis, ys);
  }


  void
  TimeModel::clearInterpolationTables()
  {
    tables.clear();
  }


  void
  TimeModel::interpolateParameters(const double DX, const double r, const double ppsi, vector<double>& output)
  {
    // psi in [0,180] deg + left-right-symmetry assumed in showers
    const double psi = fabs(utl::NormalizeAngleMinusPiPi(ppsi));

    if (interpMode == 0) {
      vector<double> intp_pars;
      for (unsigned int i = 0, n = tables.size(); i < n; ++i) {
        IntTable& t = tables[i];
        intp_pars.push_back(int2(t.psis, t.rs, t.ys, psi, r));
      }

      for (unsigned int ipar = 0; ipar < nParams; ++ipar)
        output.push_back(getShapeParameter(ipar, intp_pars, DX));
    } else {
      // warning: assuming all parameters are defined on the same grid (r,psi)
      IntTable& t0 = tables[0];
      const unsigned int ir = FindIndex(t0.rs, r);
      const unsigned int ip = FindIndex(t0.psis, psi);

      const vector<double> gp_rs(t0.rs.begin() + ir, t0.rs.begin() + ir + 2);
      const vector<double> gp_psis(t0.psis.begin() + ip, t0.psis.begin() + ip + 2);

      for (unsigned int ipar = 0; ipar < nParams; ++ipar) {
        vector<vector<double> > gpshapepars;

        for (unsigned int cur_ir = ir; cur_ir <= ir + 1; ++cur_ir) {

          vector<double> tmp;
          for (unsigned int cur_ip = ip; cur_ip <= ip + 1; ++cur_ip) {

            vector<double> gppars;  // parameters of the current grid point
            for (unsigned int i = 0, n = tables.size(); i < n; ++i)
              gppars.push_back(tables[i].ys[cur_ir][cur_ip]);
            tmp.push_back(getShapeParameter(ipar, gppars, DX));

          }
          gpshapepars.push_back(tmp);
        }

        output.push_back(int2(gp_psis, gp_rs, gpshapepars, psi, r));
      }
    }
  }


  double
  TimeModel::firstParticlePdf(const double t, const double npart)
  {
    return npart * pow((1 - cdf(t)), npart - 1) * pdf(t);
  }


  double
  TimeModel::firstParticlePdfSmeared(const double t, const double npart)
  {
    const double sigma = 25;
    const double dtau = 2;
    double f = 0;
    for (double tau = -3 * sigma; tau < 3 * sigma; tau += dtau)
      f += ROOT::Math::normal_pdf(tau, sigma, 0) * firstParticlePdf(t - tau, npart) * dtau;
    return f;
  }


  double
  TimeModel::getRisetime()
  {
    return invcdf(0.5) - invcdf(0.1);
  }


  double
  TimeModel::getFalltime()
  {
    return invcdf(0.9) - invcdf(0.5);
  }

}