Commit 6e24cd8f authored by Jan Luc Tavares's avatar Jan Luc Tavares

adicionados programas whatisit e NRAnalyser

parent 3a0fadb1
# Calculos necessarios para uma analise do tipo NRA
#
# Baseado no programa "Calibra RBS"
#
# OBS: Este programa nao tem funcionalidade geral e funciona especificamente para a situacao em que ele foi desenvolvido.
# Sinta-se livre para editar e adaptar.
#
# Licenca GPLv3.
from bokeh.plotting import figure, output_file, show
import scipy as sp
from scipy.optimize import curve_fit
import numpy as np
########################################################
####################### functions ######################
########################################################
def le_arq(endereco):
canal = []
contagens = []
arq = open(endereco)
for line in arq:
li=line.strip()
if (not li.startswith("#")) : #ignora linhas com comentarios
li = li.split()
canal.append(float(li[0]))
contagens.append(float(li[1]))
arq.close()
return np.array(canal), np.array(contagens)
def errf (x, A, B, C, D): #funcao erro
altura = 5000.0
interm = sp.special.erf(D*x - 991.0 - B)
y = A*altura*interm + C
# Observe que os parametros nao podem estar muito longe dos dados a serem fitados
return y
def linear(x, A, B):
return (A*x + B)
def linplot (x, A, B):
# por alguma razao o plot nao funciona com a linear normal.
i = 0
saida = []
for valor in x:
saida.append(A*x[i] + B)
return saida
def integra(dados, roi0, roi1):
'''
Funcao que soma os valores dos dados entre a regiao de interesse definida por roi0 e roi1.
Ou seja, calcula a area sob o grafico entre os valores roi0 e roi1.
'''
integral = 0
for i in range(roi[1] - roi[0]):
integral = integral + dados[(roi[0] + i)]
return integral
def expf(x, k, A, B):
return A*10*np.exp(-k*0.01*x) + B
########################################################
#################### definitions #######################
########################################################
colors = ["Blue", "DarkOliveGreen", "BlueViolet", "Brown", "CadetBlue", "Coral"]
# endereco dos arquivos de medidas. Estes arquivos estao no computador de Jan Luc. (este trecho exemplifica a observação de que esse código nao é genérico)
arquivos_padrao = ["/home/jan/2017-1/AMFI/Trabalho2NRA/files/F989.dat","/home/jan/2017-1/AMFI/Trabalho2NRA/files/F990.dat","/home/jan/2017-1/AMFI/Trabalho2NRA/files/F991.dat","/home/jan/2017-1/AMFI/Trabalho2NRA/files/F992.dat","/home/jan/2017-1/AMFI/Trabalho2NRA/files/F993.dat","/home/jan/2017-1/AMFI/Trabalho2NRA/files/F994.dat","/home/jan/2017-1/AMFI/Trabalho2NRA/files/F995.dat","/home/jan/2017-1/AMFI/Trabalho2NRA/files/F996.dat","/home/jan/2017-1/AMFI/Trabalho2NRA/files/F997.dat"]
arquivos_amostra = ['/home/jan/2017-1/AMFI/Trabalho2NRA/files/C989.dat', '/home/jan/2017-1/AMFI/Trabalho2NRA/files/C990.dat', '/home/jan/2017-1/AMFI/Trabalho2NRA/files/C991.dat', '/home/jan/2017-1/AMFI/Trabalho2NRA/files/C992.dat', '/home/jan/2017-1/AMFI/Trabalho2NRA/files/C993.dat', '/home/jan/2017-1/AMFI/Trabalho2NRA/files/C994.dat', '/home/jan/2017-1/AMFI/Trabalho2NRA/files/C995.dat', '/home/jan/2017-1/AMFI/Trabalho2NRA/files/C996.dat', '/home/jan/2017-1/AMFI/Trabalho2NRA/files/C997.dat', '/home/jan/2017-1/AMFI/Trabalho2NRA/files/C998.dat', '/home/jan/2017-1/AMFI/Trabalho2NRA/files/C999.dat', '/home/jan/2017-1/AMFI/Trabalho2NRA/files/C1000.dat', '/home/jan/2017-1/AMFI/Trabalho2NRA/files/C1001.dat', '/home/jan/2017-1/AMFI/Trabalho2NRA/files/C1002.dat', '/home/jan/2017-1/AMFI/Trabalho2NRA/files/C1003.dat', '/home/jan/2017-1/AMFI/Trabalho2NRA/files/C1004.dat', '/home/jan/2017-1/AMFI/Trabalho2NRA/files/C1005.dat', '/home/jan/2017-1/AMFI/Trabalho2NRA/files/C1006.dat']
roi = [380, 460]
#altura = 13000 #valor necessario para o ajuste da errf nao divergir
########################################################
################ Graficando o NRA ######################
########################################################
#
# output_file("grafNRA.html")
# p = figure(title="Grafico de NRA", x_axis_label='canal', y_axis_label='contagens')
# canais, contagens = le_arq(arquivos_padrao[2])
#p.line(canais, contagens, legend = "Grafico de NRA", line_width=2, line_color = colors[-1])
# show(p)
########################################################
############## Calculando as areas ####################
########################################################
# Para o padrão
i = 0
integrais_p = []
energias_p = []
for arquivo in arquivos_padrao:
canais, contagens = le_arq(arquivo)
integrais_p.append(integra(contagens, roi[0], roi[1]))
energias_p.append(int(arquivo[-7:-4]))
print (integrais_p, energias_p)
# Para a amostra
integrais = []
energias= []
for arquivo in arquivos_amostra:
canais, contagens = le_arq(arquivo)
integrais.append(integra(contagens, roi[0], roi[1]))
if int(arquivo[-7:-4]) > 900:
energias.append(int(arquivo[-7:-4]))
else:
energias.append(1000 + int(arquivo[-7:-4]))
print (integrais, energias)
########################################################
########### Graficando areas/energias #################
########################################################
output_file("areas.html")
q = figure(title="Grafico de NRA", x_axis_label='Energia do feixe', y_axis_label='Area do pico')
q.line(energias_p, integrais_p, legend = "Padrao de Al", line_width=2, line_color = colors[-1])
q.line(energias, integrais, legend = "Amostra de AlSb", line_width=2, line_color = colors[-2])
show(q)
########################################################
########################################################
######### Fitting para encontrar a borda ###############
########################################################
# Para o padrao
popt, pcov = curve_fit(errf, energias_p[0:-5], integrais_p[0:-5])
corte_p = (991 + popt[1])/popt[3]
energias_p = np.array(energias_p)
calculado_p = errf(energias_p, popt[0], popt[1], popt[2], popt[3])
# Para a amostra
popt, pcov = curve_fit(errf, energias[0:-5], integrais[0:-5])
corte = (991 + popt[1])/popt[3]
energias = np.array(energias)
calculado = errf(energias, popt[0], popt[1], popt[2], popt[3])
output_file("fittings.html")
p = figure(title="Ajuste das curvas usando erfc", x_axis_label='energia nominal do feixe incidente', y_axis_label='area do pico')
p.line(energias_p, integrais_p, legend = "Padrao de Al", line_width=2, line_color = colors[-3])
p.line(energias_p, calculado_p, legend = "Ajuste para o padrão", line_width=2, line_color = colors[-4])
p.line(energias, integrais, legend = "Amostra AlSb", line_width=2, line_color = colors[1])
p.line(energias, calculado, legend = "Ajuste para a amostra", line_width=2, line_color = colors[0])
show(p)
########################################################
######### Fitting para determinar de/dx ################
########################################################
energia, dedx = le_arq("/home/jan/2017-1/AMFI/Trabalho2NRA/H-em-Al.txt")
popt, pcov = curve_fit(expf, energia, dedx)
dedx_Al=expf(992,popt[0], popt[1], popt[2])
energia, dedx = le_arq("/home/jan/2017-1/AMFI/Trabalho2NRA/H-em-AlSb.txt")
popt, pcov = curve_fit(expf, energia, dedx)
dedx_AlSb=expf(992,popt[0], popt[1], popt[2])
########################################################
############ Conversao para profundidade ###############
########################################################
prof_p = (energias_p - corte_p)/dedx_Al
prof = (energias - corte)/dedx_AlSb
# Graficando a conversao
output_file("areas_profundidade.html")
q = figure(title="Grafico de NRA", x_axis_label='Profundidade aproximada de penetracao do feixe (Å)', y_axis_label='Area do pico')
q.line(prof_p, integrais_p, legend = "Padrao de Al", line_width=2, line_color = colors[-1])
q.line(prof, integrais, legend = "Amostra de AlSb", line_width=2, line_color = colors[-2])
show(q)
########################################################
############ Conversao para concentracao ###############
########################################################
bg=0
conc = (np.array(integrais) - bg)*(dedx_Al/dedx_AlSb)*100/(max(integrais_p))
type(integrais)
type(energias)
# Graficando a conversao
output_file("areas_profundidade_concentracaorel.html")
r = figure(title="Grafico de NRA", x_axis_label='Profundidade aproximada de penetracao do feixe (Å)', y_axis_label='Concentracao Relativa')
#q.line(prof_p, integrais_p, legend = "Padrao de Al", line_width=2, line_color = colors[-1])
r.line(prof, conc, legend = "Convertido para conc. relativa", line_width=2, line_color = colors[-2])
show(r)
0 0
1 0
2 0
3 0
4 0
5 0
6 0
7 0
8 0
9 0
10 0
11 0
12 0
13 0
14 0
15 0
16 0
17 0
18 0
19 0
20 0
21 0
22 0
23 0
24 0
25 0
26 0
27 0
28 0
29 0
30 0
31 0
32 0
33 0
34 0
35 0
36 0
37 0
38 407
39 426
40 417
41 394
42 416
43 395
44 388
45 380
46 336
47 318
48 354
49 341
50 410
51 349
52 423
53 425
54 518
55 507
56 614
57 592
58 530
59 481
60 451
61 344
62 299
63 282
64 265
65 294
66 314
67 364
68 416
69 447
70 465
71 451
72 395
73 322
74 265
75 201
76 150
77 121
78 119
79 120
80 95
81 101
82 100
83 110
84 104
85 95
86 122
87 114
88 100
89 116
90 106
91 74
92 84
93 94
94 82
95 91
96 82
97 84
98 95
99 118
100 100
101 112
102 118
103 90
104 117
105 108
106 91
107 74
108 63
109 72
110 52
111 49
112 49
113 57
114 60
115 54
116 55
117 60
118 47
119 46
120 55
121 54
122 56
123 54
124 42
125 46
126 37
127 53
128 52
129 32
130 53
131 45
132 36
133 40
134 37
135 54
136 32
137 47
138 39
139 40
140 43
141 34
142 34
143 40
144 34
145 39
146 46
147 36
148 60
149 39
150 46
151 44
152 44
153 50
154 35
155 34
156 37
157 30
158 47
159 38
160 58
161 33
162 38
163 35
164 33
165 39
166 35
167 31
168 36
169 42
170 43
171 41
172 43
173 56
174 48
175 40
176 66
177 48
178 46
179 50
180 71
181 47
182 66
183 53
184 49
185 38
186 49
187 50
188 47
189 59
190 41
191 39
192 30
193 27
194 34
195 51
196 33
197 43
198 28
199 32
200 38
201 30
202 36
203 31
204 29
205 42
206 32
207 41
208 42
209 35
210 38
211 24
212 32
213 20
214 27
215 34
216 28
217 28
218 28
219 42
220 30
221 24
222 42
223 29
224 28
225 35
226 26
227 33
228 24
229 32
230 27
231 34
232 37
233 33
234 52
235 37
236 38
237 37
238 46
239 35
240 46
241 40
242 33
243 43
244 31
245 36
246 35
247 43
248 35
249 22
250 42
251 22
252 27
253 28
254 23
255 33
256 32
257 31
258 35
259 30
260 35
261 26
262 30
263 22
264 36
265 23
266 22
267 36
268 27
269 21
270 30
271 34
272 29
273 32
274 29
275 28
276 39
277 26
278 25
279 27
280 31
281 26
282 30
283 31
284 34
285 29
286 34
287 36
288 29
289 35
290 49
291 41
292 46
293 39
294 42
295 33
296 39
297 42
298 31
299 31
300 43
301 34
302 41
303 41
304 39
305 42
306 56
307 33
308 58
309 63
310 49
311 46
312 53
313 51
314 60
315 54
316 53
317 41
318 55
319 58
320 44
321 38
322 29
323 26
324 30
325 36
326 34
327 26
328 21
329 30
330 13
331 14
332 17
333 19
334 13
335 12
336 7
337 15
338 12
339 4
340 17
341 12
342 16
343 9
344 11
345 16
346 13
347 14
348 13
349 15
350 16
351 18
352 12
353 13
354 17
355 17
356 14
357 30
358 15
359 13
360 21
361 15
362 9
363 20
364 22
365 22
366 14
367 19
368 16
369 16
370 17
371 21
372 17
373 12
374 14
375 33
376 11
377 19
378 22
379 14
380 25
381 18
382 21
383 23
384 22
385 26
386 13
387 17
388 19
389 29
390 25
391 19
392 23
393 29
394 24
395 24
396 21
397 26
398 32
399 29
400 34
401 26
402 38
403 21
404 38
405 40
406 33
407 24
408 36
409 29
410 46
411 20
412 36
413 32
414 46
415 56
416 48
417 41
418 43
419 49
420 36
421 49
422 37
423 52
424 59
425 39
426 40
427 63
428 34
429 44
430 51
431 48
432 47
433 42
434 24
435 38
436 39
437 35
438 25
439 14
440 17
441 28
442 21
443 16
444 15
445 14
446 11
447 9
448 6
449 7
450 2
451 7
452 7
453 2
454 4
455 2
456 5
457 2
458 2
459 4
460 2
461 2
462 3
463 4
464 4
465 1
466 5
467 5
468 4
469 2
470 7
471 2
472 2
473 6
474 2
475 0
476 7
477 7
478 3
479 5
480 2
481 6
482 6
483 5
484 7
485 5
486 4
487 6
488 1
489 2
490 7
491 3
492 1
493 9
494 1
495 4
496 2
497 1
498 5
499 4
500 1
501 0
502 5
503 7
504 4
505 2
506 0
507 3
508 0
509 2
510 4
511 2
0 0
1 0
2 0
3 0
4 0
5 0
6 0
7 0
8 0
9 0
10 0
11 0
12 0
13 0
14 0
15 0
16 0
17 0
18 0
19 0
20 0
21 0
22 0
23 0
24 0
25 0
26 0
27 0
28 0
29 0
30 0