ALICE_RAA.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ALICE Masterclass R$_{AA}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    Laden der erforderlichen Pakete"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ALICE_RAA_Tools as at\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import pandas as pd\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    Liste der verfügbaren Zentralitäten:\n",
    " \n",
    "     0-5%, 0-10%, 5-10%, 10-20%, 20-30%, 30-40%, 40-50%, 50-60%, 60-70%, 70-80%, 80-90%\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "listCentralities = [(60,70), (70,80)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    Erstellen eines Dictionaries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dictCentralities = at.create_dictionary(listCentralities)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Event-Analyse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 Einlesen der Events"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    Zählen der Events"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "event_lines = at.count_lines('event_information.csv')\n",
    "print(\"Number of events: '\" + str(event_lines) + \"'!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    Einlesen der Event-Information aus der CSV-Datei"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df_events = pd.read_csv(\"event_information.csv\", header=None, names=['eventMult', 'eventCent'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(df_events.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Bearbeiten der Daten"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    Hinzufügen von Spalten zum DataFrame für jedes Zentralitätsintervall"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for key in dictCentralities:\n",
    "    df_events[key] = False\n",
    "\n",
    "print(df_events.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    Füllen der Zentralitätsinformation und Bestimmung der Anzahl an Events pro Zentralitätsintervall"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dictEventsCent = {}\n",
    "for row in df_events.iterrows():\n",
    "    keys = at.find_centralities(dictCentralities, row[1]['eventCent'])\n",
    "    if keys:\n",
    "        for cent in keys:\n",
    "            df_events.loc[row[0], cent] = True\n",
    "            if cent in dictEventsCent:\n",
    "                dictEventsCent[cent] += 1\n",
    "            else:\n",
    "                dictEventsCent[cent] = 1\n",
    "                \n",
    "print(df_events.head())\n",
    "print(dictEventsCent)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3 Darstellung der Ergebnisse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    Erzeugung & Darstellung der Histogramme zur Teilchenmultiplizität"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hTPC_0_100, bins = at.create_event_histos(df_events['eventMult'])\n",
    "hTPC_60_70, bins = at.create_event_histos(df_events[df_events['60-70'] == True]['eventMult'])\n",
    "hTPC_70_80, bins = at.create_event_histos(df_events[df_events['70-80'] == True]['eventMult'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hTPC_0_100_err = np.sqrt(hTPC_0_100)\n",
    "hTPC_60_70_err = np.sqrt(hTPC_60_70)\n",
    "hTPC_70_80_err = np.sqrt(hTPC_70_80)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.errorbar(bins[:-1], hTPC_0_100, yerr=hTPC_0_100_err, fmt='g+', label='0-100%')\n",
    "plt.errorbar(bins[:-1], hTPC_60_70, yerr=hTPC_60_70_err, fmt='b+', label='60-70%')\n",
    "plt.errorbar(bins[:-1], hTPC_70_80, yerr=hTPC_70_80_err, fmt='r+', label='70-80%')\n",
    "#zu den Einstllungen unter 'fmt' ist folgende Seite nützlich\n",
    "#https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.plot.html#matplotlib.pyplot.plot\n",
    "\n",
    "plt.xlabel('Number of TPC tracks')\n",
    "plt.ylabel('Counts')\n",
    "plt.ylim(1,2E4)\n",
    "plt.yscale('log')\n",
    "\n",
    "plt.legend(loc='upper right')\n",
    "plt.savefig('TPC_tracks.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    2D-Plot: Teilchenspurmultiplizität gegen Zentralität"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist2d(df_events['eventMult'],df_events['eventCent'], bins=(200,100), range=[[0, 2000], [0, 100]], cmap=plt.cm.jet, norm=mpl.colors.LogNorm())\n",
    "\n",
    "plt.colorbar()\n",
    "plt.xlabel('Number of TPC tracks')\n",
    "plt.ylabel('Centrality')\n",
    "plt.savefig('TPC_tracks_2D.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Teilchenspuren und Impulsspektren"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    Liste der verfügbaren Zentralitäten:\n",
    " \n",
    "     0-5%, 0-10%, 5-10%, 10-20%, 20-30%, 30-40%, 40-50%, 50-60%, 60-70%, 70-80%, 80-90%"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.1 Einlesen der Pb-Pb Daten"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Number of tracks: '38042122'! ~600MB Datei (ungepackt)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dictMomCent = {}\n",
    "for key in dictCentralities:\n",
    "    dictMomCent[key] = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_tracks = pd.read_pickle(\"./track_info.pkl\", 'bz2').to_numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Sortieren nach Zentralität"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for line in df_tracks:\n",
    "    keys = at.find_centralities(dictCentralities, line[1])\n",
    "    if keys:\n",
    "        for key in keys:\n",
    "            dictMomCent[key].append(line[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Erzeugen der Histogramme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "binsPt = at.get_bins()\n",
    "binWidth = at.get_bin_width(binsPt)\n",
    "x_bin_width = np.asarray(binWidth)/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hPbPb_60_70, _ = np.histogram(dictMomCent['60-70'], binsPt)\n",
    "hPbPb_60_70_err = np.sqrt(hPbPb_60_70) / dictEventsCent['60-70']\n",
    "hPbPb_60_70 = hPbPb_60_70 / binWidth\n",
    "hPbPb_60_70 = hPbPb_60_70 / dictEventsCent['60-70']\n",
    "\n",
    "hPbPb_70_80, _ = np.histogram(dictMomCent['70-80'], binsPt)\n",
    "# Spezialfall für 70-80%, ein Bin ist Null - problematisch für spätere Division - daher 0 durch 1 ersetzen\n",
    "hPbPb_70_80 = np.where(hPbPb_70_80==0,1,hPbPb_70_80)\n",
    "hPbPb_70_80_err = np.sqrt(hPbPb_70_80) / dictEventsCent['70-80']\n",
    "hPbPb_70_80 = hPbPb_70_80 / binWidth\n",
    "hPbPb_70_80 = hPbPb_70_80 / dictEventsCent['70-80']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4 Einlesen der Proton-Proton Referenz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pp = np.genfromtxt('pp_reference.dat')\n",
    "print(pp)\n",
    "pp_err = np.divide(pp,10) #Annahme: 10% Fehler auf der pp Messung"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.5 Darstellung der Transversalimpulsspektren geladener Teilchen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.errorbar(binsPt[:-1], hPbPb_60_70, xerr=x_bin_width, yerr=hPbPb_60_70_err, fmt='g.', label='60-70%')\n",
    "plt.errorbar(binsPt[:-1], hPbPb_70_80, xerr=x_bin_width, yerr=hPbPb_70_80_err, fmt='b.', label='70-80%')\n",
    "plt.errorbar(binsPt[:-1], pp,          xerr=x_bin_width, yerr=pp_err,          fmt='y+', label='pp')\n",
    "\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "plt.xlim(0.1, 20)\n",
    "plt.ylim(1E-6, 1E4)\n",
    "\n",
    "plt.xlabel('$p_{T}$')\n",
    "plt.ylabel('Counts')\n",
    "plt.legend(loc='upper right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.6 Definition der Zentralitätsklassen"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Definition der Anzahl an Kollisionen für einzelne Zentralitäten"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dictNColl = {'0-5': 1686.87,\n",
    "             '0-10': 1502.7,\n",
    "             '5-10': 1319.89,\n",
    "             '10-20': 923.89,\n",
    "             '20-30': 558.68,\n",
    "             '30-40': 321.20,\n",
    "             '40-50': 171.67,\n",
    "             '50-60': 85.13,\n",
    "             '60-70': 38.51,\n",
    "             '70-80': 15.78,\n",
    "             '80-90': 6.32,\n",
    "             'pp':1.\n",
    "            }\n",
    "\n",
    "print(dictNColl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.7 Berechnung und Darstellung des R$_{\\text{CP}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RCP_60_70 = np.divide(hPbPb_60_70 / dictNColl['60-70'], hPbPb_70_80 / dictNColl['70-80'])                                                                "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RCP_60_70_err = at.fehlerberechnung(hPbPb_60_70, hPbPb_60_70_err, dictNColl['60-70'], hPbPb_70_80, hPbPb_70_80_err, dictNColl['70-80'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.errorbar(binsPt[:-1], RCP_60_70, xerr=x_bin_width, yerr=RCP_60_70_err, fmt='b.', label='60-70%')\n",
    "\n",
    "plt.ylim(-0.1,2)\n",
    "plt.xlabel('$p_{T}$')\n",
    "plt.ylabel(r'$R_{cp}=\\frac{1/N^{AA,X}_{evt}1/<N_{coll,X}>dN^{AA,X}/dp_t}{1/N^{PbPb,per}_{evt}1/<N_{coll,per}>dN^{PbPb,per}/dp_t}$')\n",
    "plt.legend(loc='upper right')\n",
    "plt.savefig('RCP_compared.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.8 Berechnung und Darstellung des R$_\\text{AA}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RAA_60_70 = np.divide(hPbPb_60_70 / dictNColl['60-70'], pp)\n",
    "RAA_70_80 = np.divide(hPbPb_70_80 / dictNColl['70-80'], pp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RAA_60_70_err = at.fehlerberechnung(hPbPb_60_70, hPbPb_60_70_err, dictNColl['60-70'], pp, pp_err, dictNColl['pp'])\n",
    "RAA_70_80_err = at.fehlerberechnung(hPbPb_70_80, hPbPb_70_80_err, dictNColl['70-80'], pp, pp_err, dictNColl['pp'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.errorbar(binsPt[:-1], RAA_60_70, xerr=x_bin_width, yerr=RAA_60_70_err, fmt='b.', label='60-70%')\n",
    "plt.errorbar(binsPt[:-1], RAA_70_80, xerr=x_bin_width, yerr=RAA_70_80_err, fmt='g.', label='70-80%')\n",
    "\n",
    "plt.ylim(0.01,1.2)\n",
    "plt.xlabel('$p_{T}$')\n",
    "plt.ylabel(r'$R_{AA}=\\frac{1/N^{AA}_{evt}dN^{AA}/dp_t}{<N_{coll}>1/N_{pp_{evt}dN^{pp}/dp_t}}$')\n",
    "plt.legend(loc='upper right')\n",
    "plt.savefig('RAA_compared.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}