Modelowanie komputerowe  procesu oddziaływania z materią ciężkich cząstek naładowanych

  Promieniowanie jądrowe jest niezwykle cennym źródłem informacji o strukturze i wielu własnościach materii. Aby te własności dobrze poznać, należy w jak najlepszym stopniu zrozumieć procesy przechodzenia promieniowania przez materię. Oczywiście wszelkie takie procesy zależą od typu promieniowania i jego energii oraz od własności materiału ośrodka. W naszym szczególnym przypadku rozpatrujemy strumienie pionów, protonów, deuteronów, cząstek alfa, oraz jąder węgla, czyli tzw. ciężkich cząstek naładowanych. (Za "cząstki ciężkie" uważamy takie, których masa jest dużo większa od masy elektronu.)

1. Prawo osłabienia wiązki promieniowania

    Samo oddziaływanie promieniowania jądrowego z materią realizuje się poprzez wiele procesów związanych z oddziaływaniem padających cząstek ze składnikami (elektronami  lub  jądrami atomowymi) materii w danym absorbencie. Prawdopodobieństwo zajścia takich procesów opisuje się wprowadzając pojęcie przekroju czynnego na dany proces  [cm²]. Przyjmijmy, iż kierujemy strumień n cząstek na pewien bardzo cienki absorbent o grubości dx, który w jednostce objętości, np. w 1 cm^3,  zawiera N centrów oddziaływania (np. atomów). Ważne jest, aby grubość dx była na tyle mała, aby jedne centra oddziaływań nie przesłaniały innych. Po przejściu naszej wiązki n cząstek niektóre z nich zostaną z niej usunięte (przez oddziaływania z centrami oddz.). Oznaczmy ten ubytek przez dn. Logiczne jest, iż dn jest wprost proporcjonalne do dx. Innymi słowy; im grubszy absorbent, tym więcej cząstek pochłania. Oprócz grubości, podobny wpływ na dn mają także inne wielkości charakteryzujące nasz układ, a mianowicie gęstość centrów oddziaływania N  oraz pierwotna liczba  n cząstek w wiązce. Powyższe rozważania zapiszmy ściśle przy uwzględnieniu wspomnianego na początku przekroju czynnego, który traktujemy tu jako współczynnik proporcjonalności:

Zauważmy, że dn, to liczba cząstek które ubyły z wiązki, stąd znak minus po prawej stronie. Przekrój czynny zapiszemy jako:

Zauważmy, że wyrażenie to zawiera w mianowniku liczbę centrów rozpraszania na jednostkę objętości N, a więc przekrój czynny odpowiada  efektywnej powierzchnia jednego centrum rozpraszającego. Musi to być powierzchnia bardzo mała i dlatego jednostką przekroju czynnego jest  1barn=10^-24 cm².  

 Skoncentrujmy się teraz na wyznaczeniu n w funkcji grubości warstwy x tj. na wyznaczeniu n(x). Zauważyć można, iż równanie (1) jest równaniem różniczkowym o zmiennych rozdzielonych i może być zapisane w postaci:

Po scałkowaniu obustronnym powyższego równania i oznaczeniu liczby cząstek, które padały na warstwę x przez n₀, otrzymujemy postać:

gdzie n(x)/n_o nazywamy współczynnikiem transmisji promieniowania przez warstwę o grubości x;  µ_x- liniowy współczynnik osłabienia wiązki; µ_d- masowy współczynnik osłabienia;  d [g/cm²]- tzw. gęstość powierzchniowa absorbentu otrzymana przez pomnożenie grubości x przez gęstość absorbentu  ( ). Warto dodać, iż jednostka g/cm² jest powszechnie używana w fizyce jądrowej, gdyż stanowi miarę grubości warstwy ośrodka, uwzględniającą zarazem jego gęstość.

Zauważmy, że rozważania nasze dotyczyły procesów, w których  cząstki są całkowicie usuwane z pierwotnej wiązki wskutek oddziaływania z atomami ośrodka. Dotyczy to przede wszystkim fotonów. Ciężkie cząstki naładowane przechodząc przez ośrodek materialny nie są usuwane z wiązki w pojedynczych aktach oddziaływania, ale stopniowo tracą swa energię w procesie strat jonizacyjnych i rozproszeń wielokrotnych. 

2. Straty energii na jonizacje atomów ośrodka

Cząstki promieniowania przechodząc przez materię oddziaływują poprzez swoje pole elektryczne z elektronami atomów wybijając je i powodując jonizację. Oczywiście w wyniku tego procesu padające cząstki tracą swoją energię. I ten właśnie proces nieelastycznych zderzeń z elektronami atomów stanowi dla ciężkich cząstek naładowanych podstawowy typ ich oddziaływania z materią (absorbentem). 

Związek między stratą energii dE cząstki, a jej przebiegiem dx w danym absorbencie określa tzw. wzór Bethe-Blocha; [1] wyrażający stosunek dE/dx tj. straty energii cząstki na jednostkę długości jej przebiegu w danym materiale:

We wzorze tym z jest liczbą atomową padającej cząstki, Z liczbą atomową absorbentu (liczbą protonów w jądrze atomu ośrodka), A liczbą masową absorbentu (sumaryczną liczbą protonów i neutronów w jądrze), c prędkością światła w próżni, K- stałą stanowiącą kombinację kilku stałych fizycznych i równą 307 keVcm²/g,   ß = v/c =pc/E - aktualną prędkością cząstki wyrażoną w jednostkach prędkości światła ( p=pęd; E=energia całkowita), pozostałe wielkości zaś:

gdzie M oznacza masę cząstki, m_e masę elektronu, a I jest średnim potencjałem jonizacji i wzbudzenia atomów absorbentu.

     Jak widać z postaci wzoru, strata energii na jednostkę drogi w absorbencie jest wprost proporcjonalna do kwadratu ładunku cząstki padającej ( z²), do stosunku liczby protonów do wszystkich nukleonów w jądrze atomu (Z/A), a co najważniejsze jest odwrotnie proporcjonalna do kwadratu prędkości cząstki  . Oznacza to, że w miarę przechodzenia cząstki przez absorbent straty energii będą bardzo szybko rosnąć. Cząstka będzie na kolejnych odcinkach swego toru deponować coraz to większe ilości energii w postaci zjonizowanych atomów, aż do zatrzymania się po przebyciu określonej drogi. Warto jednak dodać, iż wzór ten traci słuszność dla bardzo małych energii cząstek, gdyż nie uwzględnia występujących wtedy procesów chwytania przez cząstkę elektronów i ponownego ich oddawania.

    Podany wzór Bethe-Blocha uwzględnia też zjawiska relatywistyczne. Czynnik logarytmiczny, gdzie kwadrat prędkości jest w liczniku i gdzie występuje też tzw. czynnik Lorentza - gamma, wzór (6), rośnie wraz ze wzrostem prędkości. Prowadzi to do słabego zwiększania się strat jonizacyjnych dla największych energii cząstek, kiedy już prędkości bliskie są prędkości światła. W rezultacie kształt zależności strat jonizacyjnych od energii cząstek ma przebieg taki, jak pokazuje to Rys.2.4, skopiowany z pozycji [3] podanych referencji. Obserwujemy bardzo duże straty jonizacyjne dla małych energii, istnienie obszaru energii gdzie straty te maja wartości minimalne i słaby wzrost dla energii największych. Dla cząstek o różnych masach i ładunkach krzywe mają różny przebieg, co także wynika ze wzoru (5).

Pokazana na Rys 2.4 zależność strat jonizacyjnych od energii cząstki prowadzi do charakterystycznej zależności średniej gęstości jonizacji w funkcji drogi cząstki w absorbencie  (Rys. 2.5 - zaczerpnięty także z pozycji [3] podanych referencji.). Na początkowym odcinku toru straty są stosunkowo niewielkie, następnie rosną w pobliżu końca toru cząstki w materii. Natychmiast nasuwa się możliwość zastosowania tej zależności dla celów radioterapii. Cząstka przebiega odcinek zdrowej tkanki słabo na nią oddziałując, natomiast działając silnie na tkankę, którą chcemy naświetlać. Krzywa pokazana na Rys.2.5 nosi nazwę "krzywej Bragga".

3. Rozproszenia wielokrotne

Oprócz strat energii naładowanej cząstki przechodzącej przez materię obserwujemy także zjawisko wielokrotnych rozproszeń tejże cząstki, zilustrowane na Rys.3.. 

Rys 3. Wielkości opisujące wielokrotne rozproszenia cząstki przechodzącej przez absorbent o grubości  . Symbol "plane" odnosi się do rzutu toru cząstki na płaszczyznę wyznaczoną przez kierunek pierwotny lotu cząstki i dowolny kierunek do niego  prostopadły (zakłada się symetrię cylindryczną).

Oddziaływanie cząstki z materią prowadzi nie tylko do zmiany wartości prędkości cząstki, ale także do zmiany jej kierunku. Obserwujemy to jako zmianę kierunku oraz  przesunięcie wyjściowego toru cząstki w odniesieniu do toru cząstki w chwili wchodzenia do absorbentu. 

Im większe straty energetyczne następują, tym większe jest prawdopodobieństwo większego odchylenia toru cząstki. Dla opisania zjawisk wielokrotnych rozproszeń cząstki przechodzącej przez materię, stosuje się wzór [1]

gdzie oznacza odchylenie standardowe rozkładu normalnego, któremu z dobrym przybliżeniem podlega rozkład wielkości charakteryzujących rozproszenia wielokrotne;  L jest gruboscią materiału, a tzw. jednostką (długością) radiacyjną. Wizualizacja symulacji tego procesu jest dostępna w naszej aplikacji (na dole strony www).

4. Modelowanie komputerowe

    W poniżej zamieszczonym programie do modelowania zjawiska przejścia cząstek naładowanych przez zadany absorbent posłużyliśmy się podanymi wyżej wzorami. Użytkownik może określić wartości następujących wielkości: 

• początkową energię kinetyczną padającej cząstki, 

• typ absorbentu, 

• grubość absorbentu, 

• wielkość kroku dx.

Wielkość kroku ważna jest ze względu na dokładność obliczeń; im mniejsze dx tym większa dokładność obliczeń- oczywiście ciągnie to za sobą konsekwencje w postaci wolniej działającego programu) oraz typ cząstki. Efektem jest narysowanie 3 poniższych wykresów:

Wykres nr 1- przedstawia zależność E(x), gdzie E jest wyrażone w MeV, a x w g/cm². Rysowana krzywa obrazuje straty energii kinetycznej padającej cząstki w zależności od drogi przebytej w absorbencie. Pokazane poniżej przykładowe wykresy pochodzą z naszego programu przy ustawionych danych domyślnych:

Wykres nr 2- przedstawia zależność dE/dx (x). Jednostki zadane analogicznie jak poprzednio. Krzywa ta, jak już podano, nosi nazwę krzywej Bragga. Na jej podstawie określa się energie padających cząstek w radioterapii względem głębokości umiejscowienia nowotworu w ludzkim ciele tak, aby widoczny "pik" przypadał właśnie w miejscu chorym, a tkanki zdrowe były poddawane jak najmniejszym niepotrzebnym jonizacjom.

Wykres nr 3- przedstawia zależność dE/dx (E). Jednostki zadane analogicznie jak poprzednio. Oś wartości jest podana w skali logarytmicznej i ze znakiem minus dla lepszego zobrazowania zachodzącego procesu.

    Sam program działa w oparciu o następujący algorytm wykorzystujący wspomniane wzory:

    Językiem użytym jest JAVA, a sam program to aplet do pliku typu HTML. W części początkowej kodu następuje standardowa deklaracja używanych zmiennych. Następnie definiujemy stałe elementy programu (przyciski, pola tekstowe etc.- funkcja INIT). Dalej wprowadzamy pętle wyboru na poszczególne pola tekstowe i pola wyboru, by odpowiednie wartości przyporządkować zmiennym według zadanych kryteriów przez Użytkownika (np. wybór ołowiu jako absorbentu). Po pętlach wyboru przechodzimy do zdefiniowania funkcji PAINT rysującej kolejne stałe elementy dla danej symulacji (np. absorbent, ramki wykresów etc.), a następnie funkcji RYSUJ, która jest kluczową funkcją naszego programu. Początkowo przebiega tzw. symulacja wstępna (niewizualizowana), która ma na celu pomiar zasięgu cząstki oraz maksymalnej wartości funkcji dE/dx (x). Używamy w tym celu wzoru Bethe-Blocha wpisanego w pętlę typu FOR przebiegającą od x=0 przez kolejne dx, aż do osiągnięcia w sumie grubości całego absorbentu D. W pętli tej obliczamy aktualną energię przez odejmowanie kolejnych wartości dE w każdym kroku od początkowej energii kinetycznej. Po tych obliczeniach odpowiednie zmienne przechowują już znane nam: maksymalną wartość dE/dx oraz maksymalny zasięg cząstki w absorbencie (w przypadku całkowitego jej pochłonięcia i zatrzymania). Dopiero na podstawie tychże danych przechodzimy do właściwej pętli symulującej i rysującej tor cząstki i 3 wykresy. Przebiega ona dokładnie tak samo jak pętla wstępna z tą różnicą, że w przypadku zatrzymania się cząstki w absorbencie, zmienna x dąży do maksymalnego zasięgu, a nie przebiega przez całą grubość D. W czasie rzeczywistym symulowania ruchu cząstki rysowane są jednocześnie wszystkie 3 wspomniane wykresy. Po zakończeniu symulacji program oczekuje dalszych poleceń (zakończenia pracy lub wprowadzenia innych danych i dalszej symulacji). Plik z kodem źródłowym dostępny jest tym samym katalogu, z którego Użytkownik otworzył tenże dokument.

Instrukcja obsługi programu:

1) Wpisać żądane wartości parametrów i KAŻDORAZOWO ZATWIERDZAĆ naciskając "ENTER"

2) Uruchomić program

3) Nacisnąć przycisk WYCZYŚĆ po zakończeniu symulacji

4) Skok do punktu 1 w celu przeprowadzania dalszej symulacji

Jeśli poniżej zamiast programu pojawił się duży szary prostokąt, należy zainstalować program do odczytu apletów dla Windows o nazwie "jre-1_5_0-windows-i586". Program ten można ściągnąć ze strony http://nowe.pl/modules/mydownloads/singlefile.php?cid=24&lid=357

UWAGA! Jeśli poniższy program mimo to nie działa, uruchom go za pomocą aplikacji APPLETVIEWER.EXE zamieszczonej w pakiecie J2SDK. Instrukcja uruchamiająca: "appletviewer.exe start.html" (pliki muszą być w tym samym katalogu). Kliknij TUTAJ aby pobrać plik z naszym programem.

LITERATURA:

[1]  CERN Particle Physics Booklet (July 1998); str: 179 - 185, (http://pdg.lbl.gov/)

[2]  A. Strzałkowski "Wstęp do fizyki jądra atomowego", Warszawa (1979); str: 14 - 21

[3]  Ewa Skrzypczak, Zygmunt Szefliński "Wstęp do fizyki jądra atomowego i cząstek elementarnych" PWN Warszawa (2002 r.); str. 32