{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Examples for error propagation\n\nHear we use the same config in `particle_amplitude.py`\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "config_str = \"\"\"\ndecay:\n    A:\n       - [R1, B]\n       - [R2, C]\n       - [R3, D]\n    R1: [C, D]\n    R2: [B, D]\n    R3: [B, C]\n\nparticle:\n    $top:\n       A: { mass: 1.86, J: 0, P: -1}\n    $finals:\n       B: { mass: 0.494, J: 0, P: -1}\n       C: { mass: 0.139, J: 0, P: -1}\n       D: { mass: 0.139, J: 0, P: -1}\n    R1: [ R1_a, R1_b ]\n    R1_a: { mass: 0.7, width: 0.05, J: 1, P: -1}\n    R1_b: { mass: 0.5, width: 0.05, J: 0, P: +1}\n    R2: { mass: 0.824, width: 0.05, J: 0, P: +1}\n    R3: { mass: 0.824, width: 0.05, J: 0, P: +1}\n\n\"\"\"\n\nimport matplotlib.pyplot as plt\nimport yaml\n\nfrom tf_pwa.config_loader import ConfigLoader\nfrom tf_pwa.histogram import Hist1D\n\nconfig = ConfigLoader(yaml.full_load(config_str))\ninput_params = {\n    \"A->R1_a.BR1_a->C.D_total_0r\": 6.0,\n    \"A->R1_b.BR1_b->C.D_total_0r\": 1.0,\n    \"A->R2.CR2->B.D_total_0r\": 2.0,\n    \"A->R3.DR3->B.C_total_0r\": 1.0,\n}\nconfig.set_params(input_params)\n\ndata = config.generate_toy(1000)\nphsp = config.generate_phsp(10000)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "After we calculated the parameters error, we will have an error matrix `config.inv_he` (using the inverse hessain).\nIt is possible to save such matrix directly by `numpy.save` and to load it by `numpy.load`.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "config.get_params_error(data=[data], phsp=[phsp])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can use the following method to profamance the error propagation\n\n\\begin{align}\\sigma_{f} = \\sqrt{\\frac{\\partial f}{\\partial x_i} V_{ij} \\frac{\\partial f}{\\partial x_j }}\\end{align}\n\nby adding some calculation here. We need to use `tensorflow` functions instead of those of `math` or `numpy`.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import tensorflow as tf\n\nwith config.params_trans() as pt:\n    a2_r = pt[\"A->R2.CR2->B.D_total_0r\"]\n    a2_i = pt[\"A->R2.CR2->B.D_total_0r\"]\n    a2_x = a2_r * tf.cos(a2_i)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And then we can calculate the error we needed as\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(a2_x.numpy(), pt.get_error(a2_x).numpy())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can also calculate some more complex examples, such as the ratio in mass range (0.75, 0.85) over full phace space.\nEven further, we can get the error of error in the meaning of error propagation.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from tf_pwa.data import data_mask\n\nm_R2 = phsp.get_mass(\"(B, D)\")\ncut_cond = (m_R2 < 0.85) & (m_R2 > 0.75)\n\namp = config.get_amplitude()\n\nwith config.params_trans() as pt1:\n    with config.params_trans() as pt:\n        int_mc = tf.reduce_sum(amp(phsp))\n        cut_phsp = data_mask(phsp, cut_cond)\n        cut_int_mc = tf.reduce_sum(amp(cut_phsp))\n        ratio = cut_int_mc / int_mc\n    error = pt.get_error(ratio)\n\nprint(ratio.numpy(), \"+/-\", error.numpy())\nprint(error.numpy(), \"+/-\", pt1.get_error(error).numpy())"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.7.9"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}