Bethe Block

import numpy as np
import matplotlib.pyplot as plt
import math

MASSAELETTRONEMEV =0.511
UMATOMEV =931.4936148385
COSTANTE= 0.1535

def bb(ZTarget,ATarget,DensitaTarget,ZProiettile,beta_pro,gamma_pro,en_rip_pro):
  mrel=MASSAELETTRONEMEV/en_rip_pro
  etaquadro=(beta_pro*gamma_pro)**2
  WMax= (2*MASSAELETTRONEMEV*etaquadro)/(1+2*mrel*math.sqrt(1+etaquadro)+mrel**2)
  I=0
  if(ZTarget<13):
    I=(12*ZTarget+7)*10**(-6)
  else:
    I=(9.76*ZTarget+58.8*((ZTarget)**(-0.19)))*10**(-6)



  return COSTANTE*DensitaTarget*(float(ZTarget)/float(ATarget))*((ZProiettile/beta_pro)**2)*(math.log((2*MASSAELETTRONEMEV*etaquadro*WMax)/I**2)-2*beta_pro**2)



#dichiarazione variabili target
ZTarget=0
ATarget=0
DensitaTarget=0
NomeTarget=""

#dichiarazione variabili proiettile
ZProiettile=0
NomeProiettile=""
AProiettile=0
en_cin_pro=0
en_rip_pro=0
en_tot_pro=0
beta_pro=0
gamma_pro=0



#lettura file configurazione
file = open("config.txt","r")

s=file.readline()
NomeTarget=s[0:len(s)-1]
s=file.readline()
ZTarget=float(s[0:len(s)])
s=file.readline()
ATarget=float(s[0:len(s)])
s=file.readline()
DensitaTarget=float(s[0:len(s)])
s=file.readline()
NomeProiettile=s[0:len(s)-1]
s=file.readline()
ZProiettile=float(s[0:len(s)])
s=file.readline()
AProiettile=float(s[0:len(s)])
s=file.readline()
en_cin_pro=float(s[0:len(s)])
enin=en_cin_pro
s=file.readline()
mode=s[0:len(s)-1]

#inizializzazione parametri particella
en_rip_pro=AProiettile*UMATOMEV
en_tot_pro=en_cin_pro+en_rip_pro
gamma_pro=en_tot_pro/en_rip_pro
beta_pro=math.sqrt(1-1/(gamma_pro*gamma_pro))


fig,ax=plt.subplots()


if mode=='s':
  #MODALITA' S , viene valutata lo stopping power in funzione dell'energia
  #vengono letti l'energia finale da raggiungere e il numero di step in cui suddividere l'intervallo
  s=file.readline()
  en_fin=float(s[0:len(s)-1])
  s=file.readline()
  nstep=int(s[0:len(s)])

  de=(en_fin-en_cin_pro)/nstep


  energia=[]
  sp=[]

  #per nstep calcola lo stopping power e aggiorna lo stato della particella incrementando la sua energia cinetica di de.
  for i in range(nstep):
    energia.append(en_cin_pro)
    sp.append(bb(ZTarget,ATarget,DensitaTarget,ZProiettile,beta_pro,gamma_pro,en_rip_pro))
    en_cin_pro+=de
    en_tot_pro=en_cin_pro+en_rip_pro
    gamma_pro=en_tot_pro/en_rip_pro
    beta_pro=math.sqrt(1-1/(gamma_pro*gamma_pro))

  ax.plot(energia,sp)
  ax.set_title('Particelle alpha di {d} Mev su {s}'.format(d=enin,s=NomeTarget))
  ax.set_xlabel('Energia[MeV]')
  ax.set_ylabel('Sp[MeV/cm]')

elif mode=='e':
  #MODALITA' E, viene valutata l'energia residua in funzione della profondita'
  #leggo lo spessore del target e in quanti step dividere lo spessore
  s=file.readline()
  spessoreTarget=float(s[0:len(s)])
  s=file.readline()
  nstep=int(s[0:len(s)-1])

  dx=(spessoreTarget)/float(nstep)

  energia=[]
  spessore=[]

  #per nstep calcolo la perdita di energia nel tratto dx e aggiorno lo stato della particella
  for i in range(nstep):
      energia.append(en_cin_pro)
      spessore.append(i*dx)
      de=(dx*bb(ZTarget,ATarget,DensitaTarget,ZProiettile,beta_pro,gamma_pro,en_rip_pro))
      en_cin_pro-=de
      en_tot_pro=en_cin_pro+en_rip_pro
      gamma_pro=en_tot_pro/en_rip_pro
      beta_pro=math.sqrt(1-1/(gamma_pro*gamma_pro))
  ax.plot(spessore,energia)
  ax.set_title('Particelle alpha di {d} Mev su {s}'.format(d=enin,s=NomeTarget))
  ax.set_xlabel('Spessore[cm]')
  ax.set_ylabel('Energia[MeV]')

elif mode=='b':
  #MODALITA' B, viene valutata lo stopping power in funzione della profondita'
   #leggo lo spessore del target e in quanti step dividere lo spessore
  s=file.readline()
  spessoreTarget=float(s[0:len(s)])
  s=file.readline()
  nstep=int(s[0:len(s)-1])


  dx=(spessoreTarget)/float(nstep)

  sp=[]
  spessore=[]
  #per nstep calcolo la perdita di energia nel tratto dx e aggiorno lo stato della particella, ma memorizzo lo stopping power
  for i in range(nstep):
      spessore.append(i*dx)
      a=bb(ZTarget,ATarget,DensitaTarget,ZProiettile,beta_pro,gamma_pro,en_rip_pro)

      de=(dx*a)
      sp.append(a)
      en_cin_pro-=de
      en_tot_pro=en_cin_pro+en_rip_pro
      gamma_pro=en_tot_pro/en_rip_pro
      beta_pro=math.sqrt(1-1/(gamma_pro*gamma_pro))
  ax.plot(spessore,sp)
  ax.set_title('Particelle alpha di {d} Mev su {6s}'.format(d=enin,s=NomeTarget))
  ax.set_xlabel('Spessore[cm]')
  ax.set_ylabel('Sp[MeV/cm]')

plt.show()