Warning, /SimTracker/TrackerMaterialAnalysis/test/dEdxWeights.ipynb is written in an unsupported language. File is not indexed.
0001 {
0002 "cells": [
0003 {
0004 "cell_type": "markdown",
0005 "metadata": {},
0006 "source": [
0007 "# dEdx weights calculation"
0008 ]
0009 },
0010 {
0011 "cell_type": "markdown",
0012 "metadata": {},
0013 "source": [
0014 "## Input"
0015 ]
0016 },
0017 {
0018 "cell_type": "markdown",
0019 "metadata": {},
0020 "source": [
0021 "This notebook expects to find a txt file called \"VolumesZPosition.txt\". This file is coming from the _SimTracker/TrackerMaterialAnalysis_ package. Just for the record, in the old days we used to correct for the Materials with spaces in their name: \n",
0022 "\n",
0023 "sed -i -e 's/M_NEMA\\ FR4\\ plate/M_NEMA_FR4_plate/g' VolumesZPosition.txt\n",
0024 "\n",
0025 "but ok now we use HGC_G10-FR4."
0026 ]
0027 },
0028 {
0029 "cell_type": "code",
0030 "execution_count": null,
0031 "metadata": {},
0032 "outputs": [],
0033 "source": [
0034 "scenario = \"D88\""
0035 ]
0036 },
0037 {
0038 "cell_type": "code",
0039 "execution_count": null,
0040 "metadata": {},
0041 "outputs": [],
0042 "source": [
0043 "inputfile = \"VolumesZPosition_%s.txt\" % (scenario)"
0044 ]
0045 },
0046 {
0047 "cell_type": "markdown",
0048 "metadata": {},
0049 "source": [
0050 "## Some imports"
0051 ]
0052 },
0053 {
0054 "cell_type": "code",
0055 "execution_count": null,
0056 "metadata": {},
0057 "outputs": [],
0058 "source": [
0059 "from collections import OrderedDict\n",
0060 "import numpy as np\n",
0061 "import pandas as pd\n",
0062 "from IPython.display import display\n",
0063 "import matplotlib.pylab as plt\n",
0064 "import seaborn as sns"
0065 ]
0066 },
0067 {
0068 "cell_type": "markdown",
0069 "metadata": {},
0070 "source": [
0071 "## Materials Properties"
0072 ]
0073 },
0074 {
0075 "cell_type": "markdown",
0076 "metadata": {},
0077 "source": [
0078 "### Materials radiation lengths from Chris, Geant (in mm)"
0079 ]
0080 },
0081 {
0082 "cell_type": "code",
0083 "execution_count": null,
0084 "metadata": {},
0085 "outputs": [],
0086 "source": [
0087 "#In mm\n",
0088 "#-------- \n",
0089 "FromChrisMatXo = OrderedDict()\n",
0090 "FromChrisMatXo['Copper'] = 14.36\n",
0091 "FromChrisMatXo['H_Scintillator'] = 425.4\n",
0092 "FromChrisMatXo['HGC_Cables'] = 0.\n",
0093 "#FromChrisMatXo['M_NEMA_FR4_plate'] = 175.\n",
0094 "#Just duplicate the above in this case\n",
0095 "FromChrisMatXo['HGC_G10-FR4'] = 175.\n",
0096 "FromChrisMatXo['Silicon'] = 93.66\n",
0097 "FromChrisMatXo['Other'] = 0.\n",
0098 "FromChrisMatXo['Air'] = 300000.\n",
0099 "FromChrisMatXo['StainlessSteel'] = 17.35\n",
0100 "FromChrisMatXo['WCu'] = 5.122\n",
0101 "FromChrisMatXo['Lead'] = 5.612\n",
0102 "#-------- \n",
0103 "FromGeantMatXo = OrderedDict()\n",
0104 "FromGeantMatXo['Copper'] = 14.3559\n",
0105 "FromGeantMatXo['H_Scintillator'] = 425.393\n",
0106 "FromGeantMatXo['HGC_Cables'] = 66.722\n",
0107 "#FromGeantMatXo['M_NEMA_FR4_plate'] = 175.056\n",
0108 "FromGeantMatXo['HGC_G10-FR4'] = 175.056\n",
0109 "FromGeantMatXo['Silicon'] = 93.6762\n",
0110 "FromGeantMatXo['Other'] = 0.\n",
0111 "FromGeantMatXo['Air'] = 301522.\n",
0112 "FromGeantMatXo['StainlessSteel'] = 17.3555\n",
0113 "FromGeantMatXo['WCu'] = 5.1225\n",
0114 "FromGeantMatXo['Lead'] = 5.6118\n",
0115 "#-------- \n",
0116 "df = pd.DataFrame.from_dict([FromChrisMatXo, FromGeantMatXo])\n",
0117 "df = df.transpose()\n",
0118 "df.columns = ['FromChris', 'FromGeant']\n",
0119 "print( '\\033[1m' + ' Radiation Lengths in mm')\n",
0120 "df"
0121 ]
0122 },
0123 {
0124 "cell_type": "markdown",
0125 "metadata": {},
0126 "source": [
0127 "### Materials nuclear interaction lengths from Chris, Geant (in mm)"
0128 ]
0129 },
0130 {
0131 "cell_type": "code",
0132 "execution_count": null,
0133 "metadata": {},
0134 "outputs": [],
0135 "source": [
0136 "#In mm\n",
0137 "#-------- \n",
0138 "MatNucIntLength = OrderedDict()\n",
0139 "FromChrisMatNucIntLength = OrderedDict()\n",
0140 "FromChrisMatNucIntLength['Copper'] = 155.1\n",
0141 "FromChrisMatNucIntLength['H_Scintillator'] = 701.3\n",
0142 "FromChrisMatNucIntLength['HGC_Cables'] = 0.\n",
0143 "#FromChrisMatNucIntLength['M_NEMA_FR4_plate'] = 484.2\n",
0144 "#Just duplicate the above in this case\n",
0145 "FromChrisMatNucIntLength['HGC_G10-FR4'] = 484.2\n",
0146 "FromChrisMatNucIntLength['Silicon'] = 457.5\n",
0147 "FromChrisMatNucIntLength['Other'] = 0.\n",
0148 "FromChrisMatNucIntLength['Air'] = 700000\n",
0149 "FromChrisMatNucIntLength['StainlessSteel'] = 166\n",
0150 "FromChrisMatNucIntLength['WCu'] = 119.9\n",
0151 "FromChrisMatNucIntLength['Lead'] = 182.6\n",
0152 "#-------- \n",
0153 "FromGeantMatNucIntLength = OrderedDict()\n",
0154 "FromGeantMatNucIntLength['Copper'] = 155.88\n",
0155 "FromGeantMatNucIntLength['H_Scintillator'] = 700.034\n",
0156 "FromGeantMatNucIntLength['HGC_Cables'] = 393.71\n",
0157 "#FromGeantMatNucIntLength['M_NEMA_FR4_plate'] = 483.429\n",
0158 "FromGeantMatNucIntLength['HGC_G10-FR4'] = 483.429\n",
0159 "FromGeantMatNucIntLength['Silicon'] = 456.628\n",
0160 "FromGeantMatNucIntLength['Other'] = 0.\n",
0161 "FromGeantMatNucIntLength['Air'] = 704083\n",
0162 "FromGeantMatNucIntLength['StainlessSteel'] = 166.272\n",
0163 "FromGeantMatNucIntLength['WCu'] = 120.105\n",
0164 "FromGeantMatNucIntLength['Lead'] = 182.472\n",
0165 "#-------- \n",
0166 "df = pd.DataFrame.from_dict([FromChrisMatNucIntLength, FromGeantMatNucIntLength])\n",
0167 "df = df.transpose()\n",
0168 "df.columns = ['FromChris', 'FromGeant']\n",
0169 "print( '\\033[1m' + ' Nuclear interaction Lengths in mm')\n",
0170 "df"
0171 ]
0172 },
0173 {
0174 "cell_type": "markdown",
0175 "metadata": {},
0176 "source": [
0177 "### Materials dEdx from PDG, Chris, Geant (in MeV/mm)"
0178 ]
0179 },
0180 {
0181 "cell_type": "code",
0182 "execution_count": null,
0183 "metadata": {},
0184 "outputs": [],
0185 "source": [
0186 "#First in MeV gr^-1 cm^2 and at the end will convert to MeV/mm\n",
0187 "dEdx = OrderedDict()\n",
0188 "#In gr/cm^3 \n",
0189 "rho = OrderedDict()\n",
0190 "#--------\n",
0191 "#Some elements necessary to build our materials\n",
0192 "#Note: Whenever in PDG values (2018) are given in more than one state (e.g gas, liquid)\n",
0193 "#we go to the materials.xml file and choose the corresponding value. \n",
0194 "dEdx['Fe'] = 1.451\n",
0195 "dEdx['Mn'] = 1.428\n",
0196 "dEdx['Cr'] = 1.456\n",
0197 "dEdx['Ni'] = 1.468\n",
0198 "dEdx['C'] = 1.745\n",
0199 "dEdx['0'] = 1.801 \n",
0200 "dEdx['H'] = 4.034\n",
0201 "dEdx['Br'] = 1.380\n",
0202 "dEdx['W'] = 1.145\n",
0203 "dEdx['N'] = 1.813\n",
0204 "dEdx['Al'] = 1.615\n",
0205 "\n",
0206 "\n",
0207 "#-------- \n",
0208 "rho['Fe'] = 7.874 \n",
0209 "rho['Mn'] = 7.440 \n",
0210 "rho['Cr'] = 7.180\n",
0211 "rho['Ni'] = 8.902 \n",
0212 "rho['C'] = 2.265 \n",
0213 "rho['0'] = 1.43**-3 #Here we choose the materials.xml value.\n",
0214 "rho['H'] = 0.07080\n",
0215 "rho['Br'] = 3.103 \n",
0216 "rho['W'] = 19.30 \n",
0217 "rho['N'] = 0.8070 \n",
0218 "rho['Ar'] = 1.639**-3 #Here we choose the materials.xml value. \n",
0219 "rho['Al'] = 2.7\n",
0220 "\n",
0221 "#--------\n",
0222 "#Detector materials (From PDG and Geometry/CMSCommonData/data/materials.xml)\n",
0223 "#You need a measurement of the density for a compound which we take from the same file\n",
0224 "dEdx['Epoxy'] = 0.53539691*dEdx['C'] + 0.13179314*dEdx['H'] + 0.33280996*dEdx['0']\n",
0225 "rho['Epoxy'] = 1.3\n",
0226 "\n",
0227 "dEdx['Kapton'] = 0.59985105*dEdx['C'] + 0.080541353*dEdx['H'] + 0.31960759*dEdx['0']\n",
0228 "rho['Kapton'] = 1.11\n",
0229 "\n",
0230 "dEdx['Copper'] = 1.403\n",
0231 "rho['Copper'] = 8.960\n",
0232 "\n",
0233 "dEdx['H_Scintillator'] = 0.91512109*dEdx['C'] + 0.084878906*dEdx['H']\n",
0234 "rho['H_Scintillator'] = 1.032\n",
0235 "\n",
0236 "dEdx['Silicon'] = 1.664\n",
0237 "rho['Silicon'] = 2.329\n",
0238 "\n",
0239 "#dEdx['M_NEMA_FR4_plate'] = 0.18077359*dEdx['Silicon'] + 0.4056325*dEdx['0'] + 0.27804208*dEdx['C'] + 0.068442752*dEdx['H'] + 0.067109079*dEdx['Br']\n",
0240 "#rho['M_NEMA_FR4_plate'] = 1.025\n",
0241 "dEdx['HGC_G10-FR4'] = 0.18077359*dEdx['Silicon'] + 0.4056325*dEdx['0'] + 0.27804208*dEdx['C'] + 0.068442752*dEdx['H'] + 0.067109079*dEdx['Br']\n",
0242 "rho['HGC_G10-FR4'] = 1.70\n",
0243 "\n",
0244 "\n",
0245 "dEdx['Other'] = 0.\n",
0246 "rho['Other'] = 0.\n",
0247 "\n",
0248 "#dEdx['Air'] = 0.7494*dEdx['N'] + 0.2369*dEdx['0'] + 0.0129*dEdx['Ar'] + 0.0008*dEdx['H']\n",
0249 "dEdx['Air'] = 0.\n",
0250 "rho['Air'] = 1.214**-3\n",
0251 "\n",
0252 "dEdx['StainlessSteel'] = 0.6996*dEdx['Fe']+0.01*dEdx['Mn']+0.19*dEdx['Cr']+0.1*dEdx['Ni']+0.0004*dEdx['C'];\n",
0253 "rho['StainlessSteel'] = 8.02\n",
0254 "\n",
0255 "dEdx['WCu'] = 0.75*dEdx['W']+0.25*dEdx['Copper']\n",
0256 "rho['WCu'] = 14.979\n",
0257 "\n",
0258 "dEdx['Lead'] = 1.122 #Pb\n",
0259 "rho['Lead'] = 11.35 #Pb\n",
0260 "\n",
0261 "dEdx['HGC_Cables'] = 0.586*dEdx['Copper'] + 0.259*dEdx['C'] + 0.138*dEdx['0'] + 0.017*dEdx['H']\n",
0262 "rho['HGC_Cables'] = 2.68\n",
0263 "\n",
0264 "#Two new elements but at this moment they are Air\n",
0265 "dEdx['HGC_EEConnector'] = 0.\n",
0266 "rho['HGC_EEConnector'] = 1.214**-3\n",
0267 "\n",
0268 "dEdx['HGC_HEConnector'] = 0.\n",
0269 "rho['HGC_HEConnector'] = 1.214**-3\n",
0270 "\n",
0271 "dEdx['Polystyrene'] = 0.91512109*dEdx['C'] + 0.084878906*dEdx['H']\n",
0272 "rho['Polystyrene'] = 1.032\n",
0273 "\n",
0274 "#--------\n",
0275 "# Now to the calculation. First, we will loop through the composite materials, not \n",
0276 "# to mess with the elements. \n",
0277 "composite_materials = ['HGC_Cables','Epoxy','Kapton','StainlessSteel','M_NEMA_FR4_plate','H_Scintillator', 'Air', 'WCu','HGC_EEConnector','HGC_HEConnector','Polystyrene']\n",
0278 "for element in dEdx: \n",
0279 " if element not in composite_materials: continue\n",
0280 " dEdx[element] = (rho[element] * dEdx[element]) / 10.\n",
0281 "#And for the rest\n",
0282 "for element in dEdx: \n",
0283 " if element in composite_materials: continue\n",
0284 " dEdx[element] = (rho[element] * dEdx[element]) / 10.\n",
0285 "\n",
0286 "\n",
0287 "#--------\n",
0288 "FromChrisdEdx = OrderedDict()\n",
0289 "FromChrisdEdx['Copper'] = 1.26\n",
0290 "FromChrisdEdx['H_Scintillator'] = 0.395\n",
0291 "FromChrisdEdx['Epoxy'] = 0.\n",
0292 "FromChrisdEdx['Kapton'] = 0.\n",
0293 "FromChrisdEdx['HGC_Cables'] = 0.\n",
0294 "FromChrisdEdx['Al'] = 0.\n",
0295 "#FromChrisdEdx['M_NEMA_FR4_plate'] = 0.322\n",
0296 "#Just duplicate the above\n",
0297 "FromChrisdEdx['HGC_G10-FR4'] = 0.322\n",
0298 "FromChrisdEdx['Silicon'] = 0.388\n",
0299 "FromChrisdEdx['Other'] = 0.\n",
0300 "FromChrisdEdx['Air'] = 0.\n",
0301 "FromChrisdEdx['StainlessSteel'] = 1.14\n",
0302 "FromChrisdEdx['WCu'] = 1.81\n",
0303 "FromChrisdEdx['Lead'] = 1.27\n",
0304 "#-------- \n",
0305 "#The range below is from commands like: \n",
0306 "# Always recheck columns in txt, so that array is pointing to desired column. \n",
0307 "#array=($(cat VolumesZPosition.txt | grep Copper | awk '{print $7}'))\n",
0308 "#IFS=$'\\n'\n",
0309 "#Highest value\n",
0310 "#echo \"${array[*]}\" | sort -nr | head -n1\n",
0311 "#Lowest value\n",
0312 "#echo \"${array[*]}\" | sort -nr | tail -n1\n",
0313 "\n",
0314 "FromGeantdEdxWithGetDEDX = OrderedDict()\n",
0315 "FromGeantdEdxWithGetDEDX['Copper'] = \"1.14494 - 1.19191\"\n",
0316 "FromGeantdEdxWithGetDEDX['H_Scintillator'] = \"0.155602 - 0.156968\"\n",
0317 "FromGeantdEdxWithGetDEDX['HGC_Cables'] = \"0.334163 - 0.343928\"\n",
0318 "#FromGeantdEdxWithGetDEDX['M_NEMA_FR4_plate'] = \"0.200079 - 0.203959\"\n",
0319 "FromGeantdEdxWithGetDEDX['Epoxy'] = 0.\n",
0320 "FromGeantdEdxWithGetDEDX['Kapton'] = 0.\n",
0321 "FromGeantdEdxWithGetDEDX['Al'] = 0.\n",
0322 "FromGeantdEdxWithGetDEDX['HGC_G10-FR4'] = \"0.200079 - 0.203959\"\n",
0323 "FromGeantdEdxWithGetDEDX['Silicon'] = \"0.298729 - 0.349232\"\n",
0324 "FromGeantdEdxWithGetDEDX['Other'] = 0.\n",
0325 "FromGeantdEdxWithGetDEDX['Air'] = \"0.000162487 - 0.000200878\"\n",
0326 "FromGeantdEdxWithGetDEDX['StainlessSteel'] = \"1.17868 - 1.24005\"\n",
0327 "FromGeantdEdxWithGetDEDX['WCu'] = \"1.4676 - 1.54044\"\n",
0328 "FromGeantdEdxWithGetDEDX['Lead'] = \"1.28276 - 1.39512\"\n",
0329 "#-------- \n",
0330 "#The range below is from commands like: \n",
0331 "# Always recheck columns in txt, so that array is pointing to desired column. \n",
0332 "#array=($(cat VolumesZPosition.txt | grep Copper | awk '{print $9}'))\n",
0333 "#IFS=$'\\n'\n",
0334 "#Highest value\n",
0335 "#echo \"${array[*]}\" | sort -nr | head -n1\n",
0336 "#Lowest value\n",
0337 "#echo \"${array[*]}\" | sort -nr | tail -n1\n",
0338 "FromGeantdEdxWithComputeTotalDEDX = OrderedDict()\n",
0339 "FromGeantdEdxWithComputeTotalDEDX['Copper'] = \"1.25937 - 1.51496\"\n",
0340 "FromGeantdEdxWithComputeTotalDEDX['H_Scintillator'] = \"0.212843 - 0.230636\"\n",
0341 "FromGeantdEdxWithComputeTotalDEDX['HGC_Cables'] = \"0.458855 - 0.507352\"\n",
0342 "#FromGeantdEdxWithComputeTotalDEDX['M_NEMA_FR4_plate'] = \"0.319344 - 0.364818\"\n",
0343 "FromGeantdEdxWithComputeTotalDEDX['Epoxy'] = 0.\n",
0344 "FromGeantdEdxWithComputeTotalDEDX['Kapton'] = 0.\n",
0345 "FromGeantdEdxWithComputeTotalDEDX['Al'] = 0.\n",
0346 "FromGeantdEdxWithComputeTotalDEDX['HGC_G10-FR4'] = \"0.319344 - 0.364818\"\n",
0347 "FromGeantdEdxWithComputeTotalDEDX['Silicon'] = \"0.38849 - 0.442262\"\n",
0348 "FromGeantdEdxWithComputeTotalDEDX['Other'] = 0.\n",
0349 "FromGeantdEdxWithComputeTotalDEDX['Air'] = \"0.000221838 - 0.000303024\"\n",
0350 "FromGeantdEdxWithComputeTotalDEDX['StainlessSteel'] = \"1.18145 - 1.40497\"\n",
0351 "FromGeantdEdxWithComputeTotalDEDX['WCu'] = \"1.85741 - 2.18222\"\n",
0352 "FromGeantdEdxWithComputeTotalDEDX['Lead'] = \"1.28388 - 1.52941\"\n",
0353 "#-------- \n",
0354 "df = pd.DataFrame.from_dict([dEdx,FromChrisdEdx, FromGeantdEdxWithGetDEDX, FromGeantdEdxWithComputeTotalDEDX])\n",
0355 "df = df.transpose()[10:]\n",
0356 "df.columns = ['FromPDG','FromChris', 'FromGeantWithGetDEDX', 'FromGeantWithComputeTotalDEDX']\n",
0357 "print( '\\033[1m' + ' dEdx in MeV/mm')\n",
0358 "df"
0359 ]
0360 },
0361 {
0362 "cell_type": "markdown",
0363 "metadata": {},
0364 "source": [
0365 "### Choice for radiation lengths values -> Geant values"
0366 ]
0367 },
0368 {
0369 "cell_type": "code",
0370 "execution_count": null,
0371 "metadata": {},
0372 "outputs": [],
0373 "source": [
0374 "# Current choice -> Geant values\n",
0375 "#In mm\n",
0376 "MatXo = OrderedDict()\n",
0377 "MatXo['Copper'] = 14.3559\n",
0378 "MatXo['H_Scintillator'] = 425.393\n",
0379 "MatXo['HGC_Cables'] = 66.722\n",
0380 "MatXo['Epoxy'] = 315.901\n",
0381 "MatXo['Kapton'] = 365.309\n",
0382 "#MatXo['M_NEMA_FR4_plate'] = 175.056\n",
0383 "MatXo['HGC_G10-FR4'] = 175.056\n",
0384 "MatXo['Silicon'] = 93.6762\n",
0385 "MatXo['Other'] = 0.\n",
0386 "MatXo['Air'] = 301522.\n",
0387 "MatXo['StainlessSteel'] = 17.3555\n",
0388 "MatXo['WCu'] = 5.1225\n",
0389 "MatXo['Lead'] = 5.6118\n",
0390 "#Definitions below is the same as with Air at the moment\n",
0391 "MatXo['HGC_EEConnector'] = 301522.\n",
0392 "MatXo['HGC_HEConnector'] = 301522.\n",
0393 "#Polystyrene definition same as H_Scintillator at the moment \n",
0394 "MatXo['Polystyrene'] = 425.393"
0395 ]
0396 },
0397 {
0398 "cell_type": "code",
0399 "execution_count": null,
0400 "metadata": {},
0401 "outputs": [],
0402 "source": [
0403 "# Let's read the input file\n",
0404 "# Mat is prepoint material, Z is post point material start, so \n",
0405 "# upper edge of prepoint material. Etable,Efull is the energy loss in the \n",
0406 "# prepoint volume with GetDEDX and ComputeTotalDEDX. \n",
0407 "matZ = pd.read_csv(inputfile, sep=\" \", header=None, names=[\"Mat\", \"Z\", \"Eta\", \"R\", \"Etable\", \"Efull\"], index_col=False)"
0408 ]
0409 },
0410 {
0411 "cell_type": "code",
0412 "execution_count": null,
0413 "metadata": {},
0414 "outputs": [],
0415 "source": [
0416 "# We will make a new column with the physical thickness of the volumes in mm\n",
0417 "matZ[\"PhysThickInmm\"] = abs(matZ[\"Z\"].shift(-1) - matZ[\"Z\"])\n",
0418 "matZ[\"PhysThickInmm\"] = matZ[\"PhysThickInmm\"].shift(1)"
0419 ]
0420 },
0421 {
0422 "cell_type": "code",
0423 "execution_count": null,
0424 "metadata": {},
0425 "outputs": [],
0426 "source": [
0427 "# We will add a column that indicates the track the relevant volume belongs to. \n",
0428 "# The logic is that right before the next track \"PhysThickInmm\" column will be \n",
0429 "# very large. \n",
0430 "matZ[\"trackflag\"] = matZ.apply(lambda row: True if row[\"PhysThickInmm\"] < 2000. else False ,axis=1)\n",
0431 "matZ[\"tracknum\"] = (( matZ[\"trackflag\"] == False) & (matZ[\"trackflag\"].shift() == True)).cumsum()\n",
0432 "matZ[\"tracknum\"] = matZ[\"tracknum\"].shift(1)\n",
0433 "matZ = matZ.drop('trackflag', 1)"
0434 ]
0435 },
0436 {
0437 "cell_type": "code",
0438 "execution_count": null,
0439 "metadata": {},
0440 "outputs": [],
0441 "source": [
0442 "# Now that we have the tracknum we will not let the last Copper volume to \n",
0443 "# its huge thickness due to track change. According to TDR BH back is at 5137.7 mm \n",
0444 "# so we put 6.0 mm for that Copper. In any case this do *not* effect the \n",
0445 "# dEdx weights calculation since it is after the sensitive material of the last layer. \n",
0446 "# UPDATE: Seems that now last volume is SS. No effect but in any case we put 35/66 mm -> 35\n",
0447 "matZ.loc[ matZ[\"PhysThickInmm\"] > 2000. , \"PhysThickInmm\" ] = 35.0"
0448 ]
0449 },
0450 {
0451 "cell_type": "code",
0452 "execution_count": null,
0453 "metadata": {},
0454 "outputs": [],
0455 "source": [
0456 "# The line below is when comparing with Chris which has no HGC_Cables. \n",
0457 "# matZ = matZ.query('Mat != \"HGC_Cables\"')\n",
0458 "\n",
0459 "# Again, adding a new column with the physical thickness of the volumes in radiation lengths\n",
0460 "matZ[\"PhysThickInXo\"] = matZ.apply(lambda row: row[\"PhysThickInmm\"] / MatXo[row[\"Mat\"]],axis=1)"
0461 ]
0462 },
0463 {
0464 "cell_type": "code",
0465 "execution_count": null,
0466 "metadata": {},
0467 "outputs": [],
0468 "source": [
0469 "# Adding a new column with the dEdx of the material that the volumes is build\n",
0470 "matZ[\"dEdx\"] = matZ.apply(lambda row: dEdx[row[\"Mat\"]],axis=1)"
0471 ]
0472 },
0473 {
0474 "cell_type": "code",
0475 "execution_count": null,
0476 "metadata": {},
0477 "outputs": [],
0478 "source": [
0479 "# Another column with the dEdx times thickness to help us with the calculation of the \n",
0480 "# final dEdx weights \n",
0481 "matZ[\"dEdxtimesdx\"] = matZ[\"dEdx\"] * matZ[\"PhysThickInmm\"]"
0482 ]
0483 },
0484 {
0485 "cell_type": "code",
0486 "execution_count": null,
0487 "metadata": {},
0488 "outputs": [],
0489 "source": [
0490 "# And here a column with the cumulative sum\n",
0491 "matZ[\"dEdxtimesdxCum\"] = matZ.groupby('tracknum')[\"dEdxtimesdx\"].cumsum()"
0492 ]
0493 },
0494 {
0495 "cell_type": "code",
0496 "execution_count": null,
0497 "metadata": {},
0498 "outputs": [],
0499 "source": [
0500 "# And here a column with the cumulative sum for Etable\n",
0501 "matZ[\"EtableCum\"] = matZ.groupby('tracknum')[\"Etable\"].cumsum()"
0502 ]
0503 },
0504 {
0505 "cell_type": "code",
0506 "execution_count": null,
0507 "metadata": {},
0508 "outputs": [],
0509 "source": [
0510 "# And here a column with the cumulative sum for Efull\n",
0511 "matZ[\"EfullCum\"] = matZ.groupby('tracknum')[\"Efull\"].cumsum()"
0512 ]
0513 },
0514 {
0515 "cell_type": "code",
0516 "execution_count": null,
0517 "metadata": {},
0518 "outputs": [],
0519 "source": [
0520 "# And here a column for the cumulative sum in radiation length \n",
0521 "matZ[\"PhysThickInXoCum\"] = matZ.groupby('tracknum')[\"PhysThickInXo\"].cumsum()"
0522 ]
0523 },
0524 {
0525 "cell_type": "code",
0526 "execution_count": null,
0527 "metadata": {},
0528 "outputs": [],
0529 "source": [
0530 "# We will add a column that indicates the layer the relevant volume belongs to. \n",
0531 "# The logic is that if the previous material is Silicon or Scintillator\n",
0532 "# we change layer.\n",
0533 "matZ[\"layerflag\"] = matZ.apply(lambda row: False if row[\"Mat\"] == \"Silicon\" or row[\"Mat\"] == \"H_Scintillator\" else True ,axis=1)\n",
0534 "matZ[\"layer\"] = ( matZ[\"layerflag\"] == True) & (matZ[\"layerflag\"].shift(1) == False) \n",
0535 "matZ[\"layer\"] = matZ.groupby('tracknum')[\"layer\"].cumsum()\n",
0536 "#The convention is that layers starts from 1\n",
0537 "matZ[\"layer\"] = matZ.apply(lambda row: row[\"layer\"] + 1,axis=1)"
0538 ]
0539 },
0540 {
0541 "cell_type": "code",
0542 "execution_count": null,
0543 "metadata": {},
0544 "outputs": [],
0545 "source": [
0546 "# Drop auxillary columns\n",
0547 "#We need layerflag for not counting the silicon/scintillator in the dedx. \n",
0548 "#matZ = matZ.drop('layerflag', 1)"
0549 ]
0550 },
0551 {
0552 "cell_type": "code",
0553 "execution_count": null,
0554 "metadata": {},
0555 "outputs": [],
0556 "source": [
0557 "# This was my misunderstanding. There should be 53 in the layers column since \n",
0558 "# after the last scintillator or silicon the index will increase, there is \n",
0559 "# material there. In case we want to filter them out uncomment the following. \n",
0560 "# Be careful! Filter chooses and not disgards!\n",
0561 "# matZ = matZ.groupby('tracknum').filter(lambda g: ~(g['layer'] == 53.0).any() ) "
0562 ]
0563 },
0564 {
0565 "cell_type": "code",
0566 "execution_count": null,
0567 "metadata": {},
0568 "outputs": [],
0569 "source": [
0570 "#display(matZ)\n",
0571 "#matZ[(matZ[\"layer\"] == 1) & matZ[\"tracknum\"] == 20]\n",
0572 "pd.set_option('display.max_columns', None)\n",
0573 "pd.set_option('display.max_rows', None)\n",
0574 "#matZ[(matZ[\"tracknum\"] == 6.0)]\n",
0575 "matZ = matZ.dropna()\n",
0576 "#matZ[ (matZ['tracknum'] == 6.0) & (matZ[\"layer\"] == 46) ]\n",
0577 "#matZ[(matZ[\"Mat\"] == \"Kapton\") & (matZ[\"tracknum\"] == 2.0)]\n",
0578 "#matZ[(matZ[\"Mat\"] == \"HGC_Cables\")]\n",
0579 "#matZ[(matZ[\"Mat\"] == \"HGC_Cables\") & (matZ[\"layer\"] == 43)]\n",
0580 "#matZ[(matZ[\"Mat\"] == \"Air\") & (abs(matZ[\"eta\"]) > 1.5) & (abs(matZ[\"eta\"]) < 2.0)]\n",
0581 "#matZ[ (matZ[\"Z\"] > 3650) & (matZ[\"Z\"] < 3700) & (matZ[\"Mat\"] == \"Silicon\") & (matZ[\"R\"] < 900)]\n",
0582 "#matZ[ (matZ[\"Z\"] > 3594) & (matZ[\"Z\"] < 3700) & (matZ[\"tracknum\"] == 4)]\n",
0583 "#matZ[ (matZ[\"Mat\"] == \"Copper\") & (matZ[\"Z\"] > 3600) & (matZ[\"Z\"] < 3750) & (matZ[\"tracknum\"] == 4) ]\n",
0584 "#matZ[ (matZ[\"Mat\"] == \"Copper\") & (matZ[\"Z\"] > 4000) & (matZ[\"Z\"] < 4750) & (matZ[\"tracknum\"] == 4) ]\n",
0585 "#matZ[ (matZ[\"Z\"] > 4000) & (matZ[\"Z\"] < 4750) & (matZ[\"tracknum\"] == 7) ]\n",
0586 "#matZ[ (matZ[\"Mat\"] == \"Lead\") & (matZ[\"tracknum\"] == 4) ]\n",
0587 "#matZ[ (matZ[\"Mat\"] == \"H_Scintillator\") ]\n",
0588 "#matZ[ (matZ[\"Mat\"] == \"HGC_Cables\") ]\n",
0589 "#matZ[ (matZ[\"Z\"] > 3210) & (matZ[\"Z\"] < 3350) & (matZ[\"tracknum\"] == 4) ]\n",
0590 "#matZ[ (matZ[\"Z\"] > 3650) & (matZ[\"Z\"] < 3700) & (matZ[\"tracknum\"] == 4) ]\n",
0591 "#matZ[ (matZ[\"Z\"] > 4750) & (matZ[\"Z\"] < 5250) & (matZ[\"tracknum\"] == 18) ]\n",
0592 "#matZ[ (matZ[\"Mat\"] == \"Copper\") & (matZ[\"Z\"] > 3210) & (matZ[\"Z\"] < 3650) & (matZ[\"tracknum\"] == 4) & (matZ[\"PhysThickInmm\"] > 5. )]\n",
0593 "#matZ[ (matZ[\"Z\"] > 4550) & (matZ[\"Z\"] < 4650)& (matZ[\"tracknum\"] == 18)]\n",
0594 "#matZ[ (matZ[\"Z\"] > 4290) & (matZ[\"Z\"] < 4350) & (matZ[\"Mat\"] == \"Epoxy\")] \n",
0595 "#matZ[ (matZ[\"Z\"] > 4290) & (matZ[\"Z\"] < 4650) & (matZ[\"Mat\"] == \"H_Scintillator\")]\n",
0596 "#matZ[ (matZ[\"Z\"] > 4290) & (matZ[\"Z\"] < 4650) & (matZ[\"Mat\"] == \"Silicon\")]\n",
0597 "#matZ[ (matZ[\"Mat\"] == \"Silicon\") & (matZ[\"tracknum\"] == 6) ]\n",
0598 "#matZ[ (matZ[\"Mat\"] == \"H_Scintillator\") & (matZ[\"tracknum\"] == 18) ]\n",
0599 "#matZ[ (matZ[\"Mat\"] == \"HGC_G10-FR4\") & (matZ[\"tracknum\"] == 4) ]\n",
0600 "#matZ[ (matZ[\"Mat\"] == \"Silicon\") & (matZ[\"Z\"] < 3900) & (matZ[\"Z\"] > 0) & (matZ[\"PhysThickInmm\"] > 0.22)]\n",
0601 "#matZ[ (matZ[\"Mat\"] == \"Silicon\") & (matZ[\"layer\"] > 28) & (matZ[\"layer\"] < 36)& (matZ[\"PhysThickInmm\"] < 0.22)]\n",
0602 "#matZ[matZ[\"layer\"] >= 0]\n",
0603 "#matZ[(matZ[\"Mat\"] == \"WCu\") & (matZ[\"Z\"] > 3950) & (matZ[\"Z\"] < 4000)]"
0604 ]
0605 },
0606 {
0607 "cell_type": "code",
0608 "execution_count": null,
0609 "metadata": {},
0610 "outputs": [],
0611 "source": [
0612 "# We will make a new column with the cassete space in mm\n",
0613 "outforSS = matZ[ (matZ[\"Mat\"] == 'StainlessSteel') & (matZ[\"PhysThickInmm\"] != 2.5) & (matZ[\"tracknum\"] == 4) ]\n",
0614 "outforSS\n",
0615 "abs(outforSS[\"Z\"].shift(-1) - outforSS[\"Z\"]).shift(1)\n",
0616 "#outforSS[\"CasseteSpace\"] = abs(outforSS[\"Z\"].shift(-1) - outforSS[\"Z\"])\n",
0617 "#outforSS[\"PhysThickInmm\"] = outforSS[\"PhysThickInmm\"].shift(1)"
0618 ]
0619 },
0620 {
0621 "cell_type": "code",
0622 "execution_count": null,
0623 "metadata": {},
0624 "outputs": [],
0625 "source": [
0626 "#I will save here the thicknesses to be able to see the differences in files\n",
0627 "for themat in ['Copper','Polystyrene','HGC_EEConnector','HGC_HEConnector','H_Scintillator','HGC_Cables','Epoxy','Kapton','HGC_G10-FR4','Silicon','Other','Air','StainlessSteel','WCu','Lead']:\n",
0628 " outtofile = matZ[ (matZ[\"Mat\"] == themat) & (matZ[\"tracknum\"] == 7) ]\n",
0629 " outtofile[[\"Mat\",\"Z\",\"PhysThickInmm\",\"tracknum\",\"layer\"]].round(3).to_csv(r'/eos/user/a/apsallid/www/material_%s_%s.txt' %(themat,scenario), index=None, sep=' ', mode='w')\n",
0630 "\n",
0631 "#matZ[matZ[\"Z\"]>0]\n",
0632 "#matZ[matZ[\"PhysThickInmm\"]> 2000.]\n",
0633 "#matZ[matZ[\"tracknum\"] == 2 ]\n",
0634 "#matZ[ (matZ[\"tracknum\"] == 20) & (matZ[\"layer\"] > 10) ]\n",
0635 "#matZ[450:500]\n",
0636 "#For the cumulative material in front of layers. \n",
0637 "newmatZ = matZ[ (matZ['tracknum'] == 6.0) & (matZ['layerflag'] == True) & (matZ['layerflag'].shift(-1) == False)]\n",
0638 "newmatZ[['layer', 'PhysThickInXoCum']]\n",
0639 "newmatZ[['layer', 'PhysThickInXoCum']].round(3).to_csv(r'/eos/user/a/apsallid/www/%s.cumulative.xo'%(scenario), header=None, index=None, sep=' ', mode='w')"
0640 ]
0641 },
0642 {
0643 "cell_type": "code",
0644 "execution_count": null,
0645 "metadata": {},
0646 "outputs": [],
0647 "source": [
0648 "# The dEdx weights calculation doesn't include the sensitive material. \n",
0649 "matZ = matZ.drop(matZ[ (matZ.Mat == \"Silicon\") | (matZ.Mat == \"H_Scintillator\") ].index)"
0650 ]
0651 },
0652 {
0653 "cell_type": "code",
0654 "execution_count": null,
0655 "metadata": {},
0656 "outputs": [],
0657 "source": [
0658 "# We will add a column with the dedx weights. First dedxtimesdx sum per layer \n",
0659 "# This is with the theoretical dedx\n",
0660 "mdEdxtimesdxsumperlayer = matZ.groupby(['tracknum','layer'])[\"dEdxtimesdx\"].sum()\n",
0661 "# And now with the detailed simulation\n",
0662 "mdEdxtimesdxsumperlayer_detailedtable = matZ.groupby(['tracknum','layer'])[\"Etable\"].sum()\n",
0663 "mdEdxtimesdxsumperlayer_detailedfull = matZ.groupby(['tracknum','layer'])[\"Efull\"].sum()\n",
0664 "\n",
0665 "#type(mdEdxtimesdxsumperlayer)\n",
0666 "mdEdxtimesdxsumperlayer = mdEdxtimesdxsumperlayer.to_frame().reset_index()\n",
0667 "mdEdxtimesdxsumperlayer_detailedtable = mdEdxtimesdxsumperlayer_detailedtable.to_frame().reset_index()\n",
0668 "mdEdxtimesdxsumperlayer_detailedfull = mdEdxtimesdxsumperlayer_detailedfull.to_frame().reset_index()\n",
0669 "\n",
0670 "#mdEdxtimesdxsumperlayer\n",
0671 "#type(mdEdxtimesdxsumperlayer)"
0672 ]
0673 },
0674 {
0675 "cell_type": "code",
0676 "execution_count": null,
0677 "metadata": {},
0678 "outputs": [],
0679 "source": [
0680 "#Let's put also the accumulated energy loss\n",
0681 "mdEdxtimesdxsumperlayer[\"dEdxtimesdxCum\"] = mdEdxtimesdxsumperlayer.groupby('tracknum')[\"dEdxtimesdx\"].cumsum()\n",
0682 "mdEdxtimesdxsumperlayer_detailedtable[\"EtableCum\"] = mdEdxtimesdxsumperlayer_detailedtable.groupby('tracknum')[\"Etable\"].cumsum()\n",
0683 "mdEdxtimesdxsumperlayer_detailedfull[\"EfullCum\"] = mdEdxtimesdxsumperlayer_detailedfull.groupby('tracknum')[\"Efull\"].cumsum()"
0684 ]
0685 },
0686 {
0687 "cell_type": "code",
0688 "execution_count": null,
0689 "metadata": {},
0690 "outputs": [],
0691 "source": [
0692 "#Final weights\n",
0693 "mdEdxtimesdxsumperlayer[\"dedxweights\"] = (mdEdxtimesdxsumperlayer[\"dEdxtimesdx\"] + mdEdxtimesdxsumperlayer[\"dEdxtimesdx\"].shift(-1))/2\n",
0694 "mdEdxtimesdxsumperlayer_detailedtable[\"dedxweights_detailedsimulationtable\"] = (mdEdxtimesdxsumperlayer_detailedtable[\"Etable\"] + mdEdxtimesdxsumperlayer_detailedtable[\"Etable\"].shift(-1))/2\n",
0695 "mdEdxtimesdxsumperlayer_detailedfull[\"dedxweights_detailedsimulationfull\"] = (mdEdxtimesdxsumperlayer_detailedfull[\"Efull\"] + mdEdxtimesdxsumperlayer_detailedfull[\"Efull\"].shift(-1))/2"
0696 ]
0697 },
0698 {
0699 "cell_type": "code",
0700 "execution_count": null,
0701 "metadata": {},
0702 "outputs": [],
0703 "source": [
0704 "#Hack for the last layer\n",
0705 "mdEdxtimesdxsumperlayer.loc[mdEdxtimesdxsumperlayer[\"layer\"] == 47.0, \"dedxweights\"] = 83.328143\n",
0706 "mdEdxtimesdxsumperlayer_detailedtable.loc[mdEdxtimesdxsumperlayer_detailedtable[\"layer\"] == 52.0, \"dedxweights_detailedsimulationtable\"] = 91.800143\n",
0707 "mdEdxtimesdxsumperlayer_detailedfull.loc[mdEdxtimesdxsumperlayer_detailedfull[\"layer\"] == 52.0, \"dedxweights_detailedsimulationfull\"] = 98.971647\n",
0708 "\n",
0709 "#Drop duplicates not needed columns\n",
0710 "mdEdxtimesdxsumperlayer_detailedtable = mdEdxtimesdxsumperlayer_detailedtable.drop(['tracknum', 'layer'], axis=1)\n",
0711 "mdEdxtimesdxsumperlayer_detailedfull = mdEdxtimesdxsumperlayer_detailedfull.drop(['tracknum', 'layer'], axis=1)\n",
0712 "\n",
0713 "#mdEdxtimesdxsumperlayer[ mdEdxtimesdxsumperlayer['tracknum'] == 6.0 ]\n",
0714 "#mdEdxtimesdxsumperlayer_detailed[ mdEdxtimesdxsumperlayer_detailed['tracknum'] == 6.0 ]\n",
0715 "mdEdxtimesdxsumperlayer = pd.concat([mdEdxtimesdxsumperlayer, mdEdxtimesdxsumperlayer_detailedtable, mdEdxtimesdxsumperlayer_detailedfull], axis=1)\n",
0716 "mdEdxtimesdxsumperlayer = mdEdxtimesdxsumperlayer.dropna()\n",
0717 "\n",
0718 "#Adding two columns dedxtable/dedxtheory and dedxfull/dedxtheory\n",
0719 "mdEdxtimesdxsumperlayer[\"EtableoverEtheory\"] = mdEdxtimesdxsumperlayer[\"Etable\"]/mdEdxtimesdxsumperlayer[\"dEdxtimesdx\"]\n",
0720 "mdEdxtimesdxsumperlayer[\"EfulloverEtheory\"] = mdEdxtimesdxsumperlayer[\"Efull\"]/mdEdxtimesdxsumperlayer[\"dEdxtimesdx\"]\n",
0721 "mdEdxtimesdxsumperlayer[ mdEdxtimesdxsumperlayer['tracknum'] == 6.0 ][:-1].dropna()"
0722 ]
0723 },
0724 {
0725 "cell_type": "code",
0726 "execution_count": null,
0727 "metadata": {},
0728 "outputs": [],
0729 "source": [
0730 "mdEdxtimesdxsumperlayer[[\"layer\",\"EtableoverEtheory\",\"EfulloverEtheory\"]][ mdEdxtimesdxsumperlayer['tracknum'] == 6.0 ][:-1] "
0731 ]
0732 },
0733 {
0734 "cell_type": "code",
0735 "execution_count": null,
0736 "metadata": {},
0737 "outputs": [],
0738 "source": [
0739 "pd.options.display.float_format = '{:,.2f}'.format\n",
0740 "\n",
0741 "mdEdxtimesdxsumperlayer[[\"layer\",\"dEdxtimesdx\",\"dedxweights\",\"dedxweights_detailedsimulationtable\",\"dedxweights_detailedsimulationfull\"]][ mdEdxtimesdxsumperlayer['tracknum'] == 6.0 ][:-1] "
0742 ]
0743 },
0744 {
0745 "cell_type": "code",
0746 "execution_count": null,
0747 "metadata": {},
0748 "outputs": [],
0749 "source": [
0750 "forthededxtable = mdEdxtimesdxsumperlayer[[\"layer\",\"dEdxtimesdx\",\"dEdxtimesdxCum\",\"EtableCum\",\"EfullCum\"]][ mdEdxtimesdxsumperlayer['tracknum'] == 6.0 ][:-1] \n",
0751 "\n",
0752 "pd.options.display.float_format = '{:,.2f}'.format\n",
0753 "\n",
0754 "mdEdxtimesdxsumperlayer[[\"layer\",\"dedxweights\"]][ mdEdxtimesdxsumperlayer['tracknum'] == 6.0 ][:-1] "
0755 ]
0756 },
0757 {
0758 "cell_type": "code",
0759 "execution_count": null,
0760 "metadata": {},
0761 "outputs": [],
0762 "source": [
0763 "forthededxtable[[\"layer\",\"dEdxtimesdxCum\"]]\n",
0764 "forthededxtable[[\"layer\",\"dEdxtimesdxCum\"]].round(3).to_csv(r'/eos/user/a/apsallid/www/%s.cumulative.dedx'%(scenario), header=None, index=None, sep=' ', mode='w')\n",
0765 "forthededxweights = mdEdxtimesdxsumperlayer[[\"layer\",\"dedxweights\"]][ mdEdxtimesdxsumperlayer['tracknum'] == 6.0 ][:-1] \n",
0766 "forthededxweights[[\"layer\",\"dedxweights\"]].round(2).to_csv(r'/eos/user/a/apsallid/www/%s.dedx.weights'%(scenario), header=None, index=None, sep=' ', mode='w')\n",
0767 "forthededxperlayer = mdEdxtimesdxsumperlayer[[\"layer\",\"dEdxtimesdx\"]][ mdEdxtimesdxsumperlayer['tracknum'] == 6.0 ][:-1]\n",
0768 "forthededxperlayer[[\"layer\",\"dEdxtimesdx\"]].round(2).to_csv(r'/eos/user/a/apsallid/www/%s.dedx.perlayer'%(scenario), header=None, index=None, sep=' ', mode='w')"
0769 ]
0770 },
0771 {
0772 "cell_type": "code",
0773 "execution_count": null,
0774 "metadata": {},
0775 "outputs": [],
0776 "source": []
0777 }
0778 ],
0779 "metadata": {
0780 "kernelspec": {
0781 "display_name": "Python 3",
0782 "language": "python",
0783 "name": "python3"
0784 },
0785 "language_info": {
0786 "codemirror_mode": {
0787 "name": "ipython",
0788 "version": 3
0789 },
0790 "file_extension": ".py",
0791 "mimetype": "text/x-python",
0792 "name": "python",
0793 "nbconvert_exporter": "python",
0794 "pygments_lexer": "ipython3",
0795 "version": "3.9.6"
0796 }
0797 },
0798 "nbformat": 4,
0799 "nbformat_minor": 2
0800 }
