EWino/results/getEWinoData.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read SLHA, SModelS and micrOmegas output and store the data in a pandas DataFrame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\", message=\"numpy.dtype size changed\")\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import glob,imp,os,sys\n",
    "from pandas import json_normalize\n",
    "import pyslha\n",
    "sys.path.append('../../smodels')\n",
    "from smodels.particlesLoader import BSMList\n",
    "from smodels.share.models.SMparticles import SMList\n",
    "from smodels.theory.model import Model\n",
    "from getRefXSecs import getXSecFor\n",
    "from microParser import micro2Py\n",
    "pd.options.mode.chained_assignment = None #Disable copy warnings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = Model(BSMparticles=BSMList, SMparticles=SMList)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "slhaFolder = '../data/slha_scanRandom'\n",
    "smodelsFolder = '../data/smodels_scanRandom'\n",
    "microFolder = '../data/micromegas_scanRandom'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Convert Experimental Results list to a dictionary\n",
    "data = []\n",
    "removeFromDict = ['topologies outside the grid',\"missing topologies\",\n",
    "                  \"missing topologies with displaced decays\", 'missing topologies with prompt decays',\n",
    "                 \"Asymmetric Branches\",\"Outside Grid\",\"Missed Topologies\",\"Long Cascades\"]\n",
    "for f in glob.glob(smodelsFolder+'/*.py'):\n",
    "    smodelsDict = imp.load_source(f.replace('.py',''),f).smodelsOutput\n",
    "    for rmKey in removeFromDict:\n",
    "        if rmKey in smodelsDict:\n",
    "            smodelsDict.pop(rmKey)\n",
    "    if 'ExptRes' in smodelsDict:\n",
    "        for res in smodelsDict['ExptRes']:\n",
    "            if 'TxNames weights (fb)' in res:\n",
    "                res.pop('TxNames weights (fb)')  \n",
    "        expList = sorted(smodelsDict['ExptRes'], key=lambda pt: pt['r'],reverse=True)\n",
    "        expDict = dict([['result%i'%i,val] for i,val in enumerate(expList)])\n",
    "        smodelsDict['ExptRes'] = expDict\n",
    "    slhaFile = smodelsDict['OutputStatus']['input file']\n",
    "    dataDict = {'filename' : os.path.basename(slhaFile)}\n",
    "    dataDict.update(smodelsDict)\n",
    "    data.append(dataDict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "96629\n"
     ]
    }
   ],
   "source": [
    "print(len(data))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Convert data to flat DataFrame:\n",
    "smodelsDF = json_normalize(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Get SLHA data:\n",
    "slhaData = []\n",
    "for f in smodelsDF['filename']:\n",
    "    slhaFile = os.path.join(slhaFolder,f)\n",
    "    slha = pyslha.readSLHAFile(slhaFile)\n",
    "    massDict = dict([[str(key),abs(val)] for key,val in slha.blocks['MASS'].items() if key > 25])\n",
    "    extparDict = dict([[str(key),val] for key,val in slha.blocks['EXTPAR'].items()])\n",
    "    widthDict = dict([[str(key),val.totalwidth] for key,val in slha.decays.items() if key > 25])\n",
    "    BRsDict = {}\n",
    "    for pdg,val in slha.decays.items():\n",
    "        if not abs(pdg) in [1000024,1000023]:\n",
    "            continue\n",
    "        initialState = model.getParticlesWith(pdg=pdg)[0].label            \n",
    "        BRsDict[initialState] = {}\n",
    "        for dec in val.decays:\n",
    "            if dec.br < 0.01: continue            \n",
    "            finalState = ','.join([model.getParticlesWith(pdg=pid)[0].label for pid in sorted(dec.ids)])\n",
    "            BRsDict[initialState][finalState] = dec.br\n",
    "    xsec8TeV = dict([ [str(proc.pidsfinal).replace('[','').replace(']','').replace(',','_').replace(' ',''),\n",
    "                   max([x.value for x in proc.get_xsecs(sqrts=8000)])*1000] \n",
    "                 for proc in slha.xsections.values()  if proc.get_xsecs(sqrts=8000)])\n",
    "    xsec13TeV = dict([ [str(proc.pidsfinal).replace('[','').replace(']','').replace(',','_').replace(' ',''),\n",
    "                   max([x.value for x in proc.get_xsecs(sqrts=13000)])*1000] \n",
    "                 for proc in slha.xsections.values()  if proc.get_xsecs(sqrts=13000)])    \n",
    "    slhaDict = {'filename' : f, 'mass' : massDict, 'width' : widthDict, 'extpar' : extparDict, \n",
    "                'xsec8TeV(fb)' : xsec8TeV, 'xsec13TeV(fb)' : xsec13TeV, 'BRs' : BRsDict}\n",
    "    slhaData.append(slhaDict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Convert to DataFrame\n",
    "slhaDF = json_normalize(slhaData)\n",
    "#Add total cross-sections:\n",
    "xsecs13 = [x for x in list(slhaDF) if 'xsec13TeV' in x]\n",
    "xsecs8 = [x for x in list(slhaDF) if 'xsec8TeV' in x]\n",
    "slhaDF['totalxsec13TeV(fb)'] = slhaDF[xsecs13].sum(axis=1)\n",
    "slhaDF['totalxsec8TeV(fb)'] = slhaDF[xsecs8].sum(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Merge with SModelS DataFrame\n",
    "dataDF = slhaDF.merge(smodelsDF,how='inner')\n",
    "# print('Final number of data points:',dataDF.shape[0])\n",
    "#print(dataDF.columns.values.tolist()) #Print all columns names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Get Micromegas data:\n",
    "microData = []\n",
    "for f in smodelsDF['filename']:\n",
    "    microFile = os.path.join(microFolder,f.replace('.slha','.micro'))\n",
    "    obs_dict = micro2Py(microFile)\n",
    "    obs_dict['filename'] = f\n",
    "    microData.append(obs_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Final number of data points: 96629\n"
     ]
    }
   ],
   "source": [
    "#Convert to DataFrame\n",
    "microDF = json_normalize(microData)\n",
    "#Merge with SModelS DataFrame\n",
    "dataDF2 = microDF.merge(dataDF,how='inner')\n",
    "print('Final number of data points:',dataDF2.shape[0])\n",
    "#print(dataDF2.columns.values.tolist()) #Print all columns names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Add reference cross-section ratio for C1C1,C1N2:\n",
    "pids = [[-1000024,1000024],[1000023,1000024],[-1000024,1000023]]\n",
    "ratio = []\n",
    "for index,row in dataDF2.iterrows():    \n",
    "    xsecRefSum = 0.0\n",
    "    xsecSum = 0.0\n",
    "    for pidpair in pids:\n",
    "        mass1 = row['mass.%i'%abs(pidpair[0])] \n",
    "        mass2 = row['mass.%i'%abs(pidpair[1])] \n",
    "        mass = (mass1+mass2)/2.\n",
    "        xsecRef = max(0,1000*getXSecFor(pidpair[0],pidpair[1],mass,13.0,'wino'))\n",
    "        xsecRefSum += xsecRef\n",
    "            \n",
    "        xsecLabel = 'xsec13TeV(fb).%i_%i' %(pidpair[0],pidpair[1])\n",
    "        try:\n",
    "            xsec = row[xsecLabel]\n",
    "            if np.isnan(xsec):\n",
    "                xsec = 0.0\n",
    "        except:\n",
    "            xsec = 0.0\n",
    "        xsecSum += xsec\n",
    "    if xsecRefSum > 0.0:\n",
    "        ratio.append(xsecSum/xsecRefSum)    \n",
    "    else:\n",
    "        ratio.append(-1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataDF2['xsecRatio'] = ratio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Save DataFrame to pickle file:\n",
    "\n",
    "dataDF2.to_pickle('EWino_scanRandom.pcl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "96629\n"
     ]
    }
   ],
   "source": [
    "print(len(dataDF2))"
   ]
  },
  {
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}