ICRC2019/Statistics.hh

Go to the documentation of this file.

void poisson_uncertainty(double x, double& plow, double& pup)
{

  plow = x - TMath::ChisquareQuantile(0.16, 2 * x) / 2;
  pup  = TMath::ChisquareQuantile(0.84, 2 * (x + 1)) / 2 - x;


}

TMatrixD jacobian(TVectorD(*f)(TVectorD, string),
                  TVectorD x, TVectorD dx,
                   string stringtype)
{

  int nx = x.GetNoElements();
  TVectorD y = f(x, stringtype);
  int ny = y.GetNoElements();

  TMatrixD jacobi(ny, nx);
  jacobi.Zero();

  TVectorD e(nx);
  e.Zero();

  for (int ix = 0; ix < nx; ++ix) {
    e.Zero();
    e[ix] = 1 * dx[ix];

    TVectorD t1 = f(x + e, stringtype);
    TVectorD t2 = f(x, stringtype);
    TVectorD tmp = (f(x + e, stringtype) - f(x, stringtype));

    TVectorD der = tmp;

    for (int iy = 0; iy < ny; ++iy) {
      jacobi[iy][ix] = der[iy] / dx[ix]; // if ny > 1 else der
    }


  }

  return jacobi;

}

TMatrixD propagate_covariance(TVectorD(*f)(TVectorD,string),
                              TVectorD x, TMatrixD xcov,
                                string stringtype)
{
  int ncol = x.GetNoElements();
  TVectorD dx(ncol);

  for (int i = 0; i < ncol; ++i)
    dx[i] = xcov[i][i] ? sqrt(xcov[i][i]) * 1e-3 : 1;

  TMatrixD jacobi = jacobian(f, x, dx, stringtype);
  TMatrixD jacT(TMatrixD::kTransposed, jacobi);
  TMatrixD final = jacobi * (xcov * jacT);

  return final; //jacobi * (xcov * jacT);

}