{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# -*- coding: utf-8 -*-\n", "\"\"\"\n", "Created on Mon Apr 26 23:41:04 2021\n", "\n", "@author: Stepan\n", "\"\"\"\n", "\n", "import numpy as np\n", "#import pandas as pd\n", "import matplotlib.pyplot as plt\n", "from urllib.error import HTTPError # recognise the error stemming from missing data\n", "#import urllib\n", "import urllib.request\n", "#from urllib.request import urlopen\n", "#############################################################################################################\n", "\n", "#Define an exception which will be raised if the data is missing and stop the notebook execution\n", "class StopExecution(Exception):\n", " def _render_traceback_(self):\n", " pass\n", " \n", "shot_no = 36602 #test discharge for which the notebook will definitely work\n", "shot = shot_no\n", "identifier='D08-W0044_shot_'+str(shot)+'_450V'\n", "\n", "def get_file(shot, identifier):\n", " #Pick the discharge to analyse \n", " URL = 'http://golem.fjfi.cvut.cz/shots/{shot}/Diagnostics/TimepixDetector/TimePix-b/{identifier}.t3pa'\n", " \n", " url = URL.format(shot=shot, identifier=identifier)\n", " print(url)\n", " try:\n", " file_name_t3pa=url\n", " with urllib.request.urlopen(file_name_t3pa) as ft3pa:\n", " line = ft3pa.readline()\n", " line = line.decode('utf‐8')\n", " ft3pa.close\n", " except HTTPError:\n", " print('File not found at %s. Aborting notebook execution.' % url)\n", " raise StopExecution\n", " \n", " \n", " return file_name_t3pa\n", "\n", "def get_file_calib(name_calib):\n", " #Pick the discharge to analyse \n", " URL = 'http://golem.fjfi.cvut.cz/shots/{shot}/Diagnostics/TimepixDetector/TimePix-b/{name_calib}.txt'\n", " \n", " url = URL.format(shot=shot, name_calib=name_calib)\n", " print(url)\n", " try:\n", " file_calib=url\n", " with urllib.request.urlopen(file_calib) as calib:\n", " line = calib.readline()\n", " line = line.decode('utf‐8')\n", " calib.close\n", " except HTTPError:\n", " print('File not found at %s. Aborting notebook execution.' % url)\n", " raise StopExecution\n", " \n", " \n", " return file_calib\n", "\n", "'''\n", "def create_t3pa_cls(shot, identifier, file_t3pa_cls,T_cl, E_cl, index_cl, matrix_index_cl, ToA_cl, ToT_cl, FToA_cl, RowNo_cl, ClmNo_cl, overflow_cl,pocet_cl, Tfirst, Tlast, Etot, dT):\n", " fileName=str(identifier)+'.t3pa_cls'\n", " t3pa_cls = open(fileName,\"w\", encoding=\"utf-8\")\n", " url = urlopen('http://golem.fjfi.cvut.cz/shots/'+str(shot)+'/Diagnostics/TimepixDetector/')\n", " \n", " #with open(file_t3pa_cls, \"w\", encoding=\"utf-8\") as t3pa_cls:\n", " t3pa_cls.write('%\\n')\n", " t3pa_cls.write('% Index\tMatrix Index\t[ RowNo, ClmNo ]\tToA\tFToA\t( ToA_in_ns )\tToT\t( ToT_in_keV )\tOverflow\\n')\n", " t3pa_cls.write('\\n')\n", " for i in range(0,len(T_cl)):\n", " t3pa_cls.write('# '+str(i+1)+', Nunmasked = '+str(pocet_cl[i])+', Nmasked = 0, Ntot = '+str(pocet_cl[i])+'\\n')\n", " t3pa_cls.write('# Tfirst = '+str(Tfirst[i])+' ns, Tlast = '+str(Tlast[i])+' ns, dT = '+str(dT[i])+' ns, Etot = '+str(Etot[i])+' keV\\n')\n", " for j in range(0,pocet_cl[i]):\n", " t3pa_cls.write(str(index_cl[i][j])+'\t'+str(matrix_index_cl[i][j])+'\t[ '+str(RowNo_cl[i][j])+',\t'+str(ClmNo_cl[i][j])+' ]\t'+str(ToA_cl[i][j])+'\t'+str(FToA_cl[i][j])+'\t( '+str(T_cl[i][j])+' ns )\t'+str(ToT_cl[i][j])+'\t( '+str(E_cl[i][j])+' keV )\t'+str(overflow_cl[i][j])+'\\n')\n", " t3pa_cls.write('\\n') \n", " \n", " t3pa_cls.close\n", " \n", " for line in url:\n", " file.write(line + '\\n')\n", " file.close()\n", " return\n", "'''\n", "#################################################################################################################\n", "# -*- coding: utf-8 -*-\n", "\"\"\"\n", "Created on Sat Apr 24 14:03:01 2021\n", "\n", "@author: Stepan\n", "\"\"\"\n", "\n", "import numpy as np\n", "#import pandas as pd\n", "import matplotlib.pyplot as plt\n", "\n", "def load_calib(file_calib):\n", " with urllib.request.urlopen(file_calib) as fc:\n", " #with open(file_calib, \"r\", encoding=\"utf-8\") as fc:\n", " \n", " #### vytvoreni 2D pole\n", " calib=[] #vytvoreni 1D pole\n", " for i in range(0,256): #tj. rozsah 0-255\n", " temp = [] # docasne pole\n", " for j in range(0,256):\n", " temp.append(0) #naplneni docasneho pole 0\n", " calib.append(temp) #naplneni pole a[] docasnym polem temp\n", " #### mam pole calib 255x255 (0-255) vyplnene nulami\n", " \n", " for i in range(0,256): #nacteni calib matice do pole calib\n", " line = fc.readline()\n", " line = line.decode('utf‐8')\n", " word=line.strip().split(' ')\n", " #print(word)\n", " for j in range(0,256):\n", " calib[i][j]=float(word[j]) #i = radek, j = sloupec0\n", " #print(calib[i][j])\n", " '''\n", " if calib[i][j]==0:\n", " print(calib[i][j])\n", " print(i)\n", " print(j)\n", " '''\n", " fc.close \n", " print('load calib') \n", " return calib\n", "\n", "def load_t3pa_file(file_t3pa):\n", " index=[]\n", " matrix_index=[]\n", " ToA=[]\n", " ToT=[]\n", " FToA=[]\n", " overflow=[]\n", " #RowNo=[]\n", " #ClmNo=[]\n", " pocet_udalosti = 0\n", " with urllib.request.urlopen(file_t3pa) as ft3pa:\n", " line = ft3pa.readline()\n", " line = line.decode('utf‐8')\n", " while True:\n", " line = ft3pa.readline()\n", " line = line.decode('utf‐8')\n", " word=line.strip().split('\\t') #v t3pa souboru je oddelovac \\t\n", " #print(word)\n", " #poc=0\n", " if line == '':\n", " break\n", " \n", " \n", " index.append(word[0])\n", " matrix_index.append(word[1])\n", " ToA.append(float(word[2]))\n", " ToT.append(float(word[3]))\n", " FToA.append(float(word[4]))\n", " overflow.append(float(word[5]))\n", " pocet_udalosti = pocet_udalosti + 1\n", " #print(ToT[pocet_udalosti-1])\n", " #RowNo.append(int(int(word[1]))//int(256))\n", " #ClmNo.append(int(int(word[1]))%int(256))\n", " \n", " #print(noise_matrix_index)\n", " \n", " \n", " #word=line.strip().split(' ')\n", " ft3pa.close\n", " print('pocet udalosti: '+str(pocet_udalosti))\n", " #print(len(index))\n", " return index, matrix_index, ToA, ToT, FToA, overflow, pocet_udalosti\n", "\n", "\n", "def noise(index, matrix_index, ToA, ToT, FToA, overflow, pocet_udalosti): #tuto fci nemus9m explicitn2 volat - volam ji v fci load_t3pa\n", " pocet=int(0) #pocet sumicich pixelu\n", " konst=int(len(index)/1000)+1\n", " #print(konst)\n", " noise_matrix_index=[]\n", " for i in range(0,konst): \n", " pom = [] # pomocne pole\n", " k=0 #pomocna promenna - udava, kolik je v czklu ve skutecnosti udalosti - aby nebyla chyba 'list index out of range'\n", " for j in range(0,1001):\n", " if i*1000+j>=len(index):\n", " break\n", " pom.append(matrix_index[i*1000+j])\n", " k=k+1\n", " for m in range(0,k):\n", " count=int(0) #pocet vvyskytu stejneho matrix index behem 1000 udalosti\n", " index_=int(-1) #budu testovat, jestli pixel na ktery koukam je sumici (abych ho nezapocital 2x)\n", " \n", " for p in range(0,pocet):\n", " #index=int(p)\n", " if pom[m]==noise_matrix_index[p]:\n", " index_=p #pixel na ktery jsem uz koukal a byl sumici\n", " break\n", " \n", " if index_ >=0 and pom[m]==noise_matrix_index[index_]:\n", " continue\n", " \n", " for l in range(0,k):\n", " if pom[m]==pom[l]:\n", " count=count+1\n", " ####podminka na sumici pixely\n", " if count>=5: #kdyz se pixel vyskytne behem tisice udalosti vicekrat nez toto cislo, je sumici \n", " noise_matrix_index.append(pom[m])\n", " #noise_matrix_index[pocet]=pom[i]\n", " pocet=pocet+1\n", " pom.clear() \n", " \n", " #print('noise pixels:')\n", " #print(noise_matrix_index)\n", " #print(len(noise_matrix_index))\n", " #print(pocet) \n", " #print('pocet udalosti:'+str(len(index)))\n", " pocet_udalosti=len(index)\n", " \n", " for n in range (0,pocet_udalosti):\n", " for o in range(0,len(noise_matrix_index)):\n", " if n >=pocet_udalosti:\n", " break\n", " if(matrix_index[n]==noise_matrix_index[o]):\n", " del matrix_index[n]\n", " del index[n]\n", " del ToA[n]\n", " del ToT[n]\n", " del FToA[n]\n", " del overflow[n]\n", " #print(len(index))\n", " pocet_udalosti=pocet_udalosti-1\n", " continue\n", " \n", " #print(pocet_udalosti) \n", " #print(len(index))\n", " return pocet_udalosti,index, matrix_index, ToA, ToT, FToA, overflow\n", "\n", "def t3pa_data(pocet_udalosti,index, matrix_index, ToA, ToT, FToA, overflow):\n", " #rovnou vyhodim sumici pixely\n", " #pocet, pocet_udalosti, noise_matrix_index, index, matrix_index, ToA, ToT, FToA, overflow=noise(index, matrix_index, ToA, ToT, FToA, overflow, pocet_udalosti)\n", " pocet_udalosti,index, matrix_index, ToA, ToT, FToA, overflow=noise(index, matrix_index, ToA, ToT, FToA, overflow, pocet_udalosti)\n", " print('pocet udalosti bez sum pixelu: '+str(pocet_udalosti))\n", " RowNo=[]\n", " ClmNo=[]\n", " for i in range(0,len(matrix_index)):\n", " RowNo.append(int(int(matrix_index[i]))//int(256))\n", " ClmNo.append(int(int(matrix_index[i]))%int(256))\n", " #print('pocet udalosti bez sum pixelu: '+str(pocet_udalosti))\n", " return index, matrix_index, ToA, ToT, FToA, overflow, pocet_udalosti, RowNo, ClmNo\n", "\n", "\n", "def energy(a, b, c, t, ToT, pocet_udalosti, RowNo, ClmNo):\n", " E=[] #energy in keV\n", " print(pocet_udalosti)\n", " for i in range (0,pocet_udalosti):\n", " sqrt=float(0.0)\n", " e1=float(0.0)\n", " e2=float(0.0)\n", " #e=float(0.0)\n", " # promenna sqrt je vnitrek odmocniny\n", " sqrt=(((float(b[RowNo[i]][ClmNo[i]])-float(a[RowNo[i]][ClmNo[i]])*float(t[RowNo[i]][ClmNo[i]])-float(ToT[i]))**2) - float(4)*float(a[RowNo[i]][ClmNo[i]])*(float(b[RowNo[i]][ClmNo[i]])*float(t[RowNo[i]][ClmNo[i]])-float(c[RowNo[i]][ClmNo[i]])-float(t[RowNo[i]][ClmNo[i]])*float(ToT[i])))\n", " if float(sqrt)<=float(0):\n", " #e=float(0)\n", " E.append(float(0))\n", " else:\n", " '''\n", " V kalibracni matici a se obcas vyskytne 0 -> ve vypoctu energie\n", " je tim padem deleni nulou -> energie diverguje. Jak to vyresit?\n", " zatim polozim energii = 0 (kdyz a=0), pak se uvidi\n", " \n", " nakonec udelam limitu vyrazu energie pro a->0 (L'hopital)\n", " '''\n", " if a[RowNo[i]][ClmNo[i]]==0:\n", " e1=((float(t[RowNo[i]][ClmNo[i]]))/float(2)) + (((float(b[RowNo[i]][ClmNo[i]])-float(a[RowNo[i]][ClmNo[i]])*float(t[RowNo[i]][ClmNo[i]])-float(ToT[i]))*(-(float(t[RowNo[i]][ClmNo[i]])))) - 2*(float(b[RowNo[i]][ClmNo[i]])*float(t[RowNo[i]][ClmNo[i]])-float(c[RowNo[i]][ClmNo[i]])-float(t[RowNo[i]][ClmNo[i]])*float(ToT[i])))/(float(2)*np.sqrt(float(sqrt)))\n", " e2=((float(t[RowNo[i]][ClmNo[i]]))/float(2)) - (((float(b[RowNo[i]][ClmNo[i]])-float(a[RowNo[i]][ClmNo[i]])*float(t[RowNo[i]][ClmNo[i]])-float(ToT[i]))*(-(float(t[RowNo[i]][ClmNo[i]])))) - 2*(float(b[RowNo[i]][ClmNo[i]])*float(t[RowNo[i]][ClmNo[i]])-float(c[RowNo[i]][ClmNo[i]])-float(t[RowNo[i]][ClmNo[i]])*float(ToT[i])))/(float(2)*np.sqrt(sqrt))\n", " else:\n", " e1=((-(float(b[RowNo[i]][ClmNo[i]]) - float(a[RowNo[i]][ClmNo[i]])*float(t[RowNo[i]][ClmNo[i]])-float(ToT[i])))+np.sqrt(float(sqrt)))/(float(2)*float(a[RowNo[i]][ClmNo[i]]))\n", " e2=((-(float(b[RowNo[i]][ClmNo[i]]) - float(a[RowNo[i]][ClmNo[i]])*float(t[RowNo[i]][ClmNo[i]])-float(ToT[i])))-np.sqrt(float(sqrt)))/(float(2)*float(a[RowNo[i]][ClmNo[i]]))\n", " \n", " #print(str(e1)+' '+str(e2))\n", " if e1<0 and e2<0:\n", " E.append(float(0))\n", " if e1>=0 and e1>e2:\n", " E.append(float(e1))\n", " if e2>=0 and e2>e1:\n", " E.append(float(e2))\n", " #print(E[i]) \n", " return E\n", "\n", "def time(ToA, FToA, pocet_udalosti, RowNo, ClmNo):\n", " T=[] #time in ns\n", " for i in range (0,pocet_udalosti):\n", " time=float(0.0)\n", " time=(float(ToA[i])-((float(FToA[i])/float(16))))*float(25)\n", " T.append(float(time))\n", " #print (T[i])\n", " ''' \n", " print(len(T))\n", " print(T[len(T)-1])\n", " for i in range(len(T)-2,len(T)):\n", " print(T[i])\n", " '''\n", " return T\n", "\n", "def time_clustering(index, matrix_index, ToA, ToT, FToA, overflow, pocet_udalosti, RowNo, ClmNo, E, T):\n", " #### vytvoreni 2D pole\n", " #calib=[] #vytvoreni 1D pole\n", " dT=float(50) #ns - 2dT je casove okno ze ktereho sbiram udalosti a davam je do jednoho clusteru \n", " '''\n", " indexNew=[]\n", " matrix_indexNew=[]\n", " ToANew=[]\n", " ToTNew=[]\n", " FToANew=[]\n", " overflowNew=[]\n", " RowNoNew=[]\n", " ClmNoNew=[]\n", " ENew=[]\n", " TNew=[]\n", " '''\n", " index_tcl=[] #casove clustery\n", " matrix_index_tcl=[]\n", " ToA_tcl=[]\n", " ToT_tcl=[]\n", " FToA_tcl=[]\n", " overflow_tcl=[]\n", " RowNo_tcl=[]\n", " ClmNo_tcl=[]\n", " E_tcl=[]\n", " T_tcl=[]\n", " \n", " pocet=[] #pocet prvku v danem casovem clusteru \n", " print('Number of measurement events (without noise pixels): '+str(len(T)))\n", " #print(T[len(T)-1])\n", " for i in range(0,len(T)):\n", " \n", " ENew=[] # docasne pole\n", " indexNew=[]\n", " matrix_indexNew=[]\n", " ToANew=[]\n", " ToTNew=[]\n", " FToANew=[]\n", " overflowNew=[]\n", " RowNoNew=[]\n", " ClmNoNew=[]\n", " TNew=[]\n", " '''\n", " if T[i]==-1:\n", " continue\n", " '''\n", " \n", " if i>=len(T):\n", " continue\n", " \n", " TNew.append(float(T[i]))\n", " ENew.append(float(E[i]))\n", " indexNew.append(int(index[i]))\n", " matrix_indexNew.append(int(matrix_index[i]))\n", " ToANew.append(float(ToA[i]))\n", " ToTNew.append(float(ToT[i]))\n", " FToANew.append(float(FToA[i]))\n", " overflowNew.append(float(overflow[i]))\n", " RowNoNew.append(int(RowNo[i]))\n", " ClmNoNew.append(int(ClmNo[i]))\n", " del T[i]\n", " del E[i]\n", " del index[i]\n", " del matrix_index[i]\n", " del ToA[i]\n", " del ToT[i]\n", " del FToA[i]\n", " del overflow[i]\n", " del RowNo[i]\n", " del ClmNo[i]\n", " \n", " \n", " #T[i]=-1 # je to blby, ale zatim nevim, jak to udelat jinak - kdyz danou udalost prepisu na -1, je jiste ze \n", " # uz ji nezapocitam znovu, smazani udalosti mi zatim nefunguje a hazi chyby\n", "\n", " for k in range(0,len(TNew)):\n", " for j in range(1,len(T)):\n", " \n", " if j>=len(T):\n", " continue\n", " \n", " if T[j]>=(TNew[k]-dT) and T[j]<=(TNew[k]+dT):\n", " TNew.append(float(T[j]))\n", " ENew.append(float(E[j]))\n", " indexNew.append(int(index[j]))\n", " matrix_indexNew.append(int(matrix_index[j]))\n", " ToANew.append(float(ToA[j]))\n", " ToTNew.append(float(ToT[j]))\n", " FToANew.append(float(FToA[j]))\n", " overflowNew.append(float(overflow[j]))\n", " RowNoNew.append(int(RowNo[j]))\n", " ClmNoNew.append(int(ClmNo[j]))\n", " del T[j]\n", " del E[j]\n", " del index[j]\n", " del matrix_index[j]\n", " del ToA[j]\n", " del ToT[j]\n", " del FToA[j]\n", " del overflow[j]\n", " del RowNo[j]\n", " del ClmNo[j]\n", " continue \n", " \n", " T_tcl.append(TNew)\n", " E_tcl.append(ENew)\n", " index_tcl.append(indexNew)\n", " matrix_index_tcl.append(matrix_indexNew)\n", " ToA_tcl.append(ToANew)\n", " ToT_tcl.append(ToTNew)\n", " FToA_tcl.append(FToANew)\n", " RowNo_tcl.append(RowNoNew)\n", " ClmNo_tcl.append(ClmNoNew)\n", " overflow_tcl.append(overflowNew)\n", " pocet.append(int(len(TNew)))\n", " \n", " \n", " #jen vypis (prochazeni casovych clusteru)\n", " '''\n", " for v in range(0,len(T_tcl)):\n", " for b in range(0, pocet[v]):\n", " print('sour ['+str(v)+'] ['+str(b)+']')\n", " #print('T '+str(T_cl[v][b])+' E '+str(E_cl[v][b]))\n", " print('T '+str(T_tcl[v][b])+' E '+str(E_tcl[v][b])+' ['+str(RowNo_tcl[v][b])+','+str(ClmNo_tcl[v][b])+']')\n", " print('')\n", " '''\n", " print('time clusters end')\n", " TNew.clear()\n", " ENew.clear()\n", " RowNoNew.clear()\n", " ClmNoNew.clear()\n", " indexNew.clear()\n", " matrix_indexNew.clear()\n", " ToANew.clear()\n", " ToTNew.clear()\n", " FToANew.clear()\n", " overflowNew.clear()\n", " return T_tcl, E_tcl, index_tcl, matrix_index_tcl,ToA_tcl, ToT_tcl, FToA_tcl, RowNo_tcl, ClmNo_tcl, overflow_tcl, pocet\n", "\n", "def clustering (T_tcl, E_tcl, index_tcl, matrix_index_tcl,ToA_tcl, ToT_tcl, FToA_tcl, RowNo_tcl, ClmNo_tcl, overflow_tcl, pocet):\n", " '''\n", " for v in range(0,len(T_tcl)):\n", " for b in range(0, pocet[v]):\n", " print('sour ['+str(v)+'] ['+str(b)+']')\n", " #print('T '+str(T_cl[v][b])+' E '+str(E_cl[v][b]))\n", " print('T '+str(T_tcl[v][b])+' E '+str(E_tcl[v][b])+' ['+str(RowNo_tcl[v][b])+','+str(ClmNo_tcl[v][b])+']')\n", " print('')\n", " ''' \n", " index_cl=[] \n", " matrix_index_cl=[]\n", " ToA_cl=[]\n", " ToT_cl=[]\n", " FToA_cl=[]\n", " overflow_cl=[]\n", " RowNo_cl=[]\n", " ClmNo_cl=[]\n", " E_cl=[]\n", " T_cl=[]\n", " \n", " Tfirst=[]\n", " Tlast=[]\n", " Etot=[]\n", " dT=[]\n", " \n", " pocet_cl=[] #pocet prvku v danem clusteru\n", " \n", " for i in range(0,len(T_tcl)):\n", " for j in range(0, pocet[i]):\n", " \n", " ENew=[] # docasne pole\n", " indexNew=[]\n", " matrix_indexNew=[]\n", " ToANew=[]\n", " ToTNew=[]\n", " FToANew=[]\n", " overflowNew=[]\n", " RowNoNew=[]\n", " ClmNoNew=[]\n", " TNew=[]\n", " \n", " if T_tcl[i][j]==-1:\n", " continue\n", " TNew.append(float(T_tcl[i][j]))\n", " ENew.append(float(E_tcl[i][j]))\n", " indexNew.append(int(index_tcl[i][j]))\n", " matrix_indexNew.append(int(matrix_index_tcl[i][j]))\n", " ToANew.append(float(ToA_tcl[i][j]))\n", " ToTNew.append(float(ToT_tcl[i][j]))\n", " FToANew.append(float(FToA_tcl[i][j]))\n", " overflowNew.append(float(overflow_tcl[i][j]))\n", " RowNoNew.append(int(RowNo_tcl[i][j]))\n", " ClmNoNew.append(int(ClmNo_tcl[i][j]))\n", " T_tcl[i][j]=-1\n", " \n", " for m in range(pocet[i]):\n", " for k in range(0,len(TNew)):\n", " #for l in range(0,pocet[i]):\n", " for l in range(0,pocet[i]):\n", " #for k in range(0,len(TNew)):\n", " \n", " if T_tcl[i][l]==-1:\n", " continue\n", " if RowNo_tcl[i][l]<=RowNoNew[k]+1 and RowNo_tcl[i][l]>=RowNoNew[k]-1 and ClmNo_tcl[i][l]<=ClmNoNew[k]+1 and ClmNo_tcl[i][l]>=ClmNoNew[k]-1:\n", " TNew.append(float(T_tcl[i][l]))\n", " ENew.append(float(E_tcl[i][l]))\n", " indexNew.append(int(index_tcl[i][l]))\n", " matrix_indexNew.append(int(matrix_index_tcl[i][l]))\n", " ToANew.append(float(ToA_tcl[i][l]))\n", " ToTNew.append(float(ToT_tcl[i][l]))\n", " FToANew.append(float(FToA_tcl[i][l]))\n", " overflowNew.append(float(overflow_tcl[i][l]))\n", " RowNoNew.append(int(RowNo_tcl[i][l]))\n", " ClmNoNew.append(int(ClmNo_tcl[i][l]))\n", " T_tcl[i][l]=-1\n", " ''' \n", " #for k in range(0,len(TNew)):\n", " for l in range(0,pocet[i]):\n", " #for l in range(0,pocet[i]):\n", " for k in range(0,len(TNew)):\n", " \n", " if T_tcl[i][l]==-1:\n", " continue\n", " if RowNo_tcl[i][l]<=RowNoNew[k]+1 and RowNo_tcl[i][l]>=RowNoNew[k]-1 and ClmNo_tcl[i][l]<=ClmNoNew[k]+1 and ClmNo_tcl[i][l]>=ClmNoNew[k]-1:\n", " TNew.append(float(T_tcl[i][l]))\n", " ENew.append(float(E_tcl[i][l]))\n", " indexNew.append(int(index_tcl[i][l]))\n", " matrix_indexNew.append(int(matrix_index_tcl[i][l]))\n", " ToANew.append(float(ToA_tcl[i][l]))\n", " ToTNew.append(float(ToT_tcl[i][l]))\n", " FToANew.append(float(FToA_tcl[i][l]))\n", " overflowNew.append(float(overflow_tcl[i][l]))\n", " RowNoNew.append(int(RowNo_tcl[i][l]))\n", " ClmNoNew.append(int(ClmNo_tcl[i][l]))\n", " T_tcl[i][l]=-1\n", " ''' \n", " T_cl.append(TNew)\n", " E_cl.append(ENew)\n", " index_cl.append(indexNew)\n", " matrix_index_cl.append(matrix_indexNew)\n", " ToA_cl.append(ToANew)\n", " ToT_cl.append(ToTNew)\n", " FToA_cl.append(FToANew)\n", " RowNo_cl.append(RowNoNew)\n", " ClmNo_cl.append(ClmNoNew)\n", " overflow_cl.append(overflowNew)\n", " pocet_cl.append(int(len(TNew))) \n", " #print(T_cl)\n", " for n in range(0,len(T_cl)):\n", " Tfirst.append(float(min(T_cl[n])))\n", " Tlast.append(float(max(T_cl[n])))\n", " Etot.append(float(sum(E_cl[n])))\n", " dT.append(float(Tlast[n])-float(Tfirst[n])) \n", " ''' \n", " for v in range(0,len(T_cl)):\n", " print('Tfirst='+str(T_first[v])+' Tlast='+str(T_last[v])+' Etot='+str(Etot[v])+' dT='+str(dT[v]))\n", " for b in range(0, pocet_cl[v]):\n", " print('cluster ['+str(v)+'] ['+str(b)+']')\n", " #print('T '+str(T_cl[v][b])+' E '+str(E_cl[v][b]))\n", " print('T '+str(T_cl[v][b])+' E '+str(E_cl[v][b])+' ['+str(RowNo_cl[v][b])+','+str(ClmNo_cl[v][b])+']')\n", " print('')\n", " ''' \n", " print('clusters end') \n", " return T_cl, E_cl, index_cl, matrix_index_cl, ToA_cl, ToT_cl, FToA_cl, RowNo_cl, ClmNo_cl, overflow_cl,pocet_cl, Tfirst, Tlast, Etot, dT\n", "\n", "def file_t3pa_cls(file_t3pa_cls,T_cl, E_cl, index_cl, matrix_index_cl, ToA_cl, ToT_cl, FToA_cl, RowNo_cl, ClmNo_cl, overflow_cl,pocet_cl, Tfirst, Tlast, Etot, dT):\n", " with open(file_t3pa_cls, \"w\", encoding=\"utf-8\") as t3pa_cls:\n", " t3pa_cls.write('%\\n')\n", " t3pa_cls.write('% Index\tMatrix Index\t[ RowNo, ClmNo ]\tToA\tFToA\t( ToA_in_ns )\tToT\t( ToT_in_keV )\tOverflow\\n')\n", " t3pa_cls.write('\\n')\n", " for i in range(0,len(T_cl)):\n", " t3pa_cls.write('# '+str(i+1)+', Nunmasked = '+str(pocet_cl[i])+', Nmasked = 0, Ntot = '+str(pocet_cl[i])+'\\n')\n", " t3pa_cls.write('# Tfirst = '+str(Tfirst[i])+' ns, Tlast = '+str(Tlast[i])+' ns, dT = '+str(dT[i])+' ns, Etot = '+str(Etot[i])+' keV\\n')\n", " for j in range(0,pocet_cl[i]):\n", " t3pa_cls.write(str(index_cl[i][j])+'\t'+str(matrix_index_cl[i][j])+'\t[ '+str(RowNo_cl[i][j])+',\t'+str(ClmNo_cl[i][j])+' ]\t'+str(ToA_cl[i][j])+'\t'+str(FToA_cl[i][j])+'\t( '+str(T_cl[i][j])+' ns )\t'+str(ToT_cl[i][j])+'\t( '+str(E_cl[i][j])+' keV )\t'+str(overflow_cl[i][j])+'\\n')\n", " t3pa_cls.write('\\n') \n", " \n", " t3pa_cls.close\n", " print('file t3pa_cls')\n", " return\n", "\n", "\n", "def energy_spectrum(Tfirst, Etot): #dela histogram - energie zaznamenana v case \n", " dt = int(100) #casove okno ze ktereho sbiram udalosti (clustery)\n", " #print(len(Tfirst))\n", " T_first=(min(Tfirst))\n", " T_last=(max(Tfirst)) #posledni z Tfirst\n", " #print(T_first)\n", " \n", " Delta_T = T_last - T_first\n", " poc = int(int(Delta_T) / int(dt)) + 1 #pocet casovych oken\n", " \n", " T=[] #cas\n", " E=[] #energie\n", " for i in range(0,poc):\n", " T.append((i*dt) + dt/2)\n", " E.append(0)\n", " \n", " for j in range(0,len(Tfirst)):\n", " time_index=0\n", " time_index=int(((Tfirst[j]-T_first)/dt))\n", " for k in range(time_index-1,time_index+2):#range(0,poc):#\n", " if k<0:\n", " k=0\n", " if k>=poc:\n", " k=poc-1\n", " if float(Tfirst[j]-T_first) < (T[k] - dt / 2):\n", " continue\n", " if float(Tfirst[j]-T_first) >= (T[k] - dt / 2) and float(Tfirst[j]-T_first) <= (T[k] + dt / 2):\n", " E[k]=float(E[k])+float(Etot[j])\n", " break\n", " if float(Tfirst[j]-T_first) > (T[k] + dt / 2):\n", " continue\n", " \n", " print('time [ns] Energy [keV]')\n", " '''\n", " E_graph=[]\n", " T_graph=[]\n", " for l in range(0, poc):\n", " if E[l]!=0:\n", " #E_graph.append(E[l])\n", " #T_graph.append(T[l])\n", " print(str(T[l])+' '+str(E[l]))\n", " '''\n", " hist(T, E, figure_E_hist)\n", " #file_hist(energy_hist, E, T)\n", " return T, E\n", "'''\n", "def file_hist(file_hist, E, T): #zapise histogram do souboru\n", " with open(file_hist, \"w\", encoding=\"utf-8\") as hist:\n", " hist.write('# time [ns] Energy [keV] \\n')\n", " for i in range(0, len(T)):\n", " #if E[i]>0:\n", " hist.write(str(T[i])+' '+str(E[i])+'\\n')\n", " hist.close\n", " return\n", "'''\n", "def hist(T, E,figure_E_hist):\n", " #plt.style.use('ggplot')\n", " plt.plot(T, E)\n", " plt.xlabel('Time [ns]')\n", " plt.ylabel('Energy [keV]')\n", " plt.xlim([0, 25000000]) #25 ms\n", " #plt.xlim(0, 209125550.0)\n", " plt.ylim(0, 10000) #10 000 keV\n", " #plt.show()\n", " plt.savefig(figure_E_hist, dpi = 1000)\n", " #data.T.plot(kind='bar')\n", " #plt.ylabel('GDP per capita')\n", " return\n", "'''\n", "caliba = 'E:/Factory_Si/caliba.txt' #H03-W0051\n", "calibb = 'E:/Factory_Si/calibb.txt' #H03-W0051\n", "calibc = 'E:/Factory_Si/calibc.txt' #H03-W0051\n", "calibt = 'E:/Factory_Si/calibt.txt' #H03-W0051\n", "#t3pa = 'E:/mereni/mereni2010/data2010/H03-W0051_Am241_No1_temp20.t3pa'\n", "t3pa = 'E:/mereni/test/H03-W0051_shot_36530_450V.t3pa'\n", "#t3pa = 'E:/mereni/test/D08-W0044_No15_450V_1s.t3pa'\n", "#t3pa = 'E:/mereni/mereni3103/I05/I05-W0037_22Na_No2_-450V_120s.t3pa'\n", "t3pa_cls='E:/mereni/test/H03-W0051_shot_36530_450V.t3pa_cls'\n", "energy_hist='E:/mereni/test/energy_hist_H03-W0051_shot_36530_450V.txt'\n", "figure_E_hist='E:/mereni/test/energy_hist_H03-W0051_shot_36530_450V.png'\n", "\n", "#pocet, noise_matrix_index=noise(t3pa)\n", "\n", "\n", "#nactu jednotlive kalibracni matice - abych to nemusel delat v kazde funkci\n", "a=load_calib(caliba)\n", "b=load_calib(calibb)\n", "c=load_calib(calibc)\n", "t=load_calib(calibt)\n", "\n", "#nactu t3pa file:\n", "index, matrix_index, ToA, ToT, FToA, overflow, pocet_udalosti = load_t3pa_file(t3pa)\n", "\n", "#vyhodim z t3pa file sumici pixely a prevedu matrix_index na cislo radku a cislo sloupce (RowNo a ClmNo)\n", "index, matrix_index, ToA, ToT, FToA, overflow, pocet_udalosti, RowNo, ClmNo = t3pa_data(pocet_udalosti, index, matrix_index, ToA, ToT, FToA, overflow)\n", "\n", "\n", "\n", "#print(pocet_udalosti)\n", "E=energy(a, b, c, t, ToT, pocet_udalosti, RowNo, ClmNo)\n", "T=time(ToA, FToA, pocet_udalosti, RowNo, ClmNo)\n", "\n", "T_tcl, E_tcl, index_tcl, matrix_index_tcl,ToA_tcl, ToT_tcl, FToA_tcl, RowNo_tcl, ClmNo_tcl, overflow_tcl, pocet=time_clustering(index, matrix_index, ToA, ToT, FToA, overflow, pocet_udalosti, RowNo, ClmNo, E, T)\n", "T_cl, E_cl, index_cl, matrix_index_cl, ToA_cl, ToT_cl, FToA_cl, RowNo_cl, ClmNo_cl, overflow_cl,pocet_cl, Tfirst, Tlast, Etot, dT = clustering(T_tcl, E_tcl, index_tcl, matrix_index_tcl,ToA_tcl, ToT_tcl, FToA_tcl, RowNo_tcl, ClmNo_tcl, overflow_tcl, pocet)\n", "\n", "file_t3pa_cls(t3pa_cls,T_cl, E_cl, index_cl, matrix_index_cl, ToA_cl, ToT_cl, FToA_cl, RowNo_cl, ClmNo_cl, overflow_cl,pocet_cl, Tfirst, Tlast, Etot, dT)\n", "\n", "energy_spectrum(Tfirst, Etot)\n", "\n", "'''\n", "\n", "##################################################################################################################\n", "\n", "\n", "\n", "\n", "#soubory, ktere ctu\n", "t3pa=get_file(shot, identifier)\n", "name_calib='caliba'\n", "caliba=get_file_calib(name_calib)\n", "name_calib='calibb'\n", "calibb=get_file_calib(name_calib)\n", "name_calib='calibc'\n", "calibc=get_file_calib(name_calib)\n", "name_calib='calibt'\n", "calibt=get_file_calib(name_calib)\n", "\n", "'''\n", "caliba = 'E:/Factory_Si/caliba.txt' #H03-W0051\n", "calibb = 'E:/Factory_Si/calibb.txt' #H03-W0051\n", "calibc = 'E:/Factory_Si/calibc.txt' #H03-W0051\n", "calibt = 'E:/Factory_Si/calibt.txt' #H03-W0051\n", "'''\n", "\n", "#vytvorene soubory:\n", "#t3pa_cls='E:/mereni/test/D08-W0044_No15_450V_1s.t3pa_cls'\n", "t3pa_cls='H03-W0051_shot_'+str(shot)+'_450V.t3pa_cls'\n", "#energy_hist='E:/mereni/test/energy_hist_D08-W0044_No15_450V_1s.txt'\n", "\n", "figure_E_hist='icon-fig'\n", "\n", "#pocet, noise_matrix_index=noise(t3pa)\n", "\n", "\n", "\n", "#nactu jednotlive kalibracni matice - abych to nemusel delat v kazde funkci\n", "a=load_calib(caliba)\n", "b=load_calib(calibb)\n", "c=load_calib(calibc)\n", "t=load_calib(calibt)\n", "\n", "#nactu a urcim jednotlive hodnoty - abych to nemusel delat v kazde funkci\n", "index, matrix_index, ToA, ToT, FToA, overflow, pocet_udalosti = load_t3pa_file(t3pa)\n", "index, matrix_index, ToA, ToT, FToA, overflow, pocet_udalosti, RowNo, ClmNo = t3pa_data(pocet_udalosti,index, matrix_index, ToA, ToT, FToA, overflow)\n", "\n", "E=energy(a, b, c, t, ToT, pocet_udalosti, RowNo, ClmNo)\n", "T=time(ToA, FToA, pocet_udalosti, RowNo, ClmNo)\n", "\n", "T_tcl, E_tcl, index_tcl, matrix_index_tcl,ToA_tcl, ToT_tcl, FToA_tcl, RowNo_tcl, ClmNo_tcl, overflow_tcl, pocet=time_clustering(index, matrix_index, ToA, ToT, FToA, overflow, pocet_udalosti, RowNo, ClmNo, E, T)\n", "T_cl, E_cl, index_cl, matrix_index_cl, ToA_cl, ToT_cl, FToA_cl, RowNo_cl, ClmNo_cl, overflow_cl,pocet_cl, Tfirst, Tlast, Etot, dT = clustering(T_tcl, E_tcl, index_tcl, matrix_index_tcl,ToA_tcl, ToT_tcl, FToA_tcl, RowNo_tcl, ClmNo_tcl, overflow_tcl, pocet)\n", "\n", "file_t3pa_cls(t3pa_cls,T_cl, E_cl, index_cl, matrix_index_cl, ToA_cl, ToT_cl, FToA_cl, RowNo_cl, ClmNo_cl, overflow_cl,pocet_cl, Tfirst, Tlast, Etot, dT)\n", "\n", "energy_spectrum(Tfirst, Etot)\n", "\n", "\n" ] }, { "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.5" } }, "nbformat": 4, "nbformat_minor": 4 }