Warsztaty z Metod Terapeutycznych: Różnice pomiędzy wersjami
(Nie pokazano 4 pośrednich wersji utworzonych przez tego samego użytkownika) | |||
Linia 258: | Linia 258: | ||
Pod adresem [http://wiki.opengatecollaboration.org/index.php/Users_Guide_V7.2:Readout_parameters_for_Radiotherapy_applications:_Actors] można znaleźć różne inne typy actora, np. Energy spectrum. Daje on możliwość zapisania widma energii w formacie .root. Aby je obejrzeć należy w terminalu otworzyć program o nazwie '''root''' i napisać magiczną komendę '''TBrowser t''' (która nie jest tak bardzo magiczna, tworzy po prostu nowy obiekt o nazwie t typu TBrowser; ale roota nie będziemy się szczegółowo uczyć; aby go opuścić trzeba napisać .quit) i otworzyć plik, który nas interesuje. | Pod adresem [http://wiki.opengatecollaboration.org/index.php/Users_Guide_V7.2:Readout_parameters_for_Radiotherapy_applications:_Actors] można znaleźć różne inne typy actora, np. Energy spectrum. Daje on możliwość zapisania widma energii w formacie .root. Aby je obejrzeć należy w terminalu otworzyć program o nazwie '''root''' i napisać magiczną komendę '''TBrowser t''' (która nie jest tak bardzo magiczna, tworzy po prostu nowy obiekt o nazwie t typu TBrowser; ale roota nie będziemy się szczegółowo uczyć; aby go opuścić trzeba napisać .quit) i otworzyć plik, który nas interesuje. | ||
+ | |||
+ | ==Zadanie 4== | ||
+ | Poniższy skrypt w programie root powinien narysować histogram pozycji z singli, których pozycja z była dodatnia. Analogicznie można by narysować histogram energii singli, których pozycja z była dodatnia. Itd. Można to podejście zastosować do policzenia ilości singli (lub hitów) spełniających jakiś warunek. | ||
+ | <source lang='c'> | ||
+ | { | ||
+ | |||
+ | |||
+ | TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("kalibracja2.root"); | ||
+ | if (!f) { | ||
+ | f = new TFile("kalibracja2.root"); | ||
+ | } | ||
+ | |||
+ | |||
+ | //TFile *f = TFile::Open("benchSPECT.root"); | ||
+ | |||
+ | TTree *Singles = (TTree*)gDirectory->Get("Singles"); | ||
+ | |||
+ | //Float_t energy; | ||
+ | Float_t gPZ, sPZ; | ||
+ | |||
+ | //Singles->SetBranchAddress("energy",&energy); | ||
+ | |||
+ | Singles->SetBranchAddress("globalPosZ",&gPZ); | ||
+ | Singles->SetBranchAddress("sourcePosZ",&sPZ); | ||
+ | |||
+ | //TH1F *ener = new TH1F("ener","energia",100,0,0.2); | ||
+ | TH1F *posZ = new TH1F("pozycjaZ", "pozycja Z", 100,0,150); | ||
+ | |||
+ | Int_t nentries = Singles->GetEntries(); | ||
+ | Int_t nbytes = 0; | ||
+ | |||
+ | |||
+ | for (Int_t i=0; i<nentries;i++) { | ||
+ | nbytes += Singles->GetEntry(i); | ||
+ | |||
+ | if (sPZ>0) | ||
+ | posZ.Fill(gPZ); | ||
+ | // ener.Fill(energy); | ||
+ | |||
+ | |||
+ | } | ||
+ | |||
+ | |||
+ | // **************************************** Plots ********************************************** | ||
+ | |||
+ | |||
+ | TCanvas cont("contours","kontury",100,100,800,600); | ||
+ | |||
+ | cont.SetFillColor(0); | ||
+ | |||
+ | posZ->SetFillColor(2); | ||
+ | posZ->SetFillStyle(3023); | ||
+ | //posZ->GetXaxis()->SetTitle(""); | ||
+ | posZ.Draw(); | ||
+ | |||
+ | } | ||
+ | |||
+ | </source> |
Aktualna wersja na dzień 13:22, 18 sty 2017
Warsztaty z Metod Terapeutycznych Rok akademicki 2016/17
Instalowanie i uruchamianie vGate 7.2
- Wejdź do katalogu "Warsztaty z metod terapeutycznych" na pulpicie. Rozpakuj plik "vGATE72.vdi.7z" (jeśli nie jest już rozpakowany).
- Uruchom "Oracle VM VirualBox".
- Dodaj nową wirtualną maszynę (przycisk "Nowa"). Nadaj jej nazwę vgate7. Wybierz typ systemu operacyjnego Linux, wersja Ubuntu (32-bit).
- W kolejnym oknie wybierz rozmiar pamięci przydzielonej dla wirtualnej maszyny 1024 MB.
- W kolejnym oknie wybierz "Użyj istniejącego pliku dysku twardego" i wskaż ścieżkę do rozpakowanego pliku "vGATE72.vdi".
- Uruchamianie wirtualnej maszyny
- Wybierz dodaną ostatnio maszynę (jeśli jest więcej niż jedna) i naciśnij "Uruchom"
- Gdy uruchomi się linux, należy wpisać hasło "virtual"
- Uruchamianie programu "Gate"
- Uruchom terminal.
- Za pomocą komend "cd nazwa_katalogu" wejdź do katalogu, w którym znajdują się Twoje skrypty. Jeśli takiego nie masz, to stwórz go w katalogu ~/Progs/Gate/gate_v7.2/
- Przydatne linki
Zadanie 1
Zdefiniuj geometrię kolimatora: płaskiego pudła ołowianego, w którym zostało wywiercone 25 okrągłych otworów z powietrza w geometrii sześciokątnej. (za 2 punkty)
Skrypty wywołuje się komendą
/control/execute nazwa_skryptu.mac
Poniżej skrypt, który definiuje właściwości świata i umieszcza w nim kostkę Rubika (3x3 pudełka utworzone z powietrza). (Na stronie [4] można znaleźć przykłady i więcej opisów jak budować geometrię)
# na początku trzeba wczytać bazę danych dostępnych materiałów
/gate/geometry/setMaterialDatabase ../GateMaterials.db
# Aby oglądać tworzoną geometrię na etapie jej definiowania;
# podczas dłuższych obliczeń należałoby z tego zrezygnować.
/vis/enable
/control/execute vis.mac
# wymiary świata
/gate/world/geometry/setXLength 100. cm
/gate/world/geometry/setYLength 100. cm
/gate/world/geometry/setZLength 50. cm
# definicja wymiarów, położenia, materiałów z których wykonano pudło
/gate/world/daughters/name pudlo
/gate/world/daughters/insert box
/gate/pudlo/geometry/setXLength 1. cm
/gate/pudlo/geometry/setYLength 1. cm
/gate/pudlo/geometry/setZLength 1. cm
/gate/pudlo/placement/setTranslation 20. 0. 0. cm
/gate/pudlo/setMaterial Air
/gate/pudlo/vis/forceWireframe
# powtórzenie pudła po 3 razy w każdą stronę poprzez dodanie do niego powtarzacza cubicArray
/gate/pudlo/repeaters/insert cubicArray
/gate/pudlo/cubicArray/setRepeatNumberX 3
/gate/pudlo/cubicArray/setRepeatNumberY 3
/gate/pudlo/cubicArray/setRepeatNumberZ 3
/gate/pudlo/cubicArray/setRepeatVector 1.1 1.1 1.1 cm
# przydatne często, gdy wprowadziło się nowe definicje geometrii, które nie zostały jeszcze zwizualizowane
/gate/geometry/rebuild
Odwołuje się przy tym do skryptu vis.mac
# V I S U A L I S A T I O N
/vis/open OGLSX # biblioteka graficzna
/vis/viewer/set/viewpointThetaPhi 60 60 # punkt obserwacji
/vis/viewer/zoom 1.5 # powiększenie
/vis/drawVolume
/vis/viewer/flush
/tracking/verbose 0 # poziom gadatliwości, czyli ile informacji wyświetlać w terminalu podczas obliczeń
/tracking/storeTrajectory 1 # ile trajektorii zapamiętywać
/vis/scene/add/trajectories
/vis/scene/endOfEventAction accumulate
Aby włożyć do pudła cylindryczny otwór
/gate/pudlo/daughters/name dziura # utworzenie otworu
/gate/pudlo/daughters/insert cylinder
/gate/dziura/geometry/setHeight 1. cm
/gate/dziura/geometry/setRmax .15 cm
/gate/dziura/geometry/setRmin 0. cm
/gate/dziura/placement/setRotationAxis 0 1 0 # cylinder można obrócić
/gate/dziura/placement/setRotationAngle 90 deg
/gate/dziura/setMaterial Air
Oczywiście w kolimatorze to nie pudło, a otwór powinien zostac odpowiednią ilość razy powtórzony.
Zadanie 2
Napromieniuj dwoma punktowymi źródełkami Technetu 99 fantom wodny poprzez kolimator i bez kolimatora. Odczytaj wyniki używając programu ImageJ. (za 2 punkty)
Ustawienie symulowanych zjawisk fizycznych może być dokonane w osobnym pliku (np. fizyka.mac) Propozycję jego zawartości można znaleźć na stronie [5] Na tej samej stronie widnieje propozycja ustawienia odcięć w ciele pacjenta ("cuts", czyli zasięg cząstek poniżej którego przestajemy śledzić cząstkę, bo ma zbyt małą energię, by była dla nas interesująca).
Do zdefiniowanego wcześniej fantomu można dołączyć "aktora", czyli zdolność interakcji z symulacją - w tym wypadku zapisania do pliku mapy rozkładu zdeponowanej energii w fantomie. Fantom możemy pociąć na mniejsze elementy (woksele), traktując go jako macierz (np. 100x100 elementów).
#=====================================================
# MATRIX FOR DOSE MAP OUTPUT
#=====================================================
/gate/actor/addActor DoseActor doseDistribution
/gate/actor/doseDistribution/save nazwa_katalogu/plik_wynikowy.hdr
/gate/actor/doseDistribution/attachTo moj_fantom
/gate/actor/doseDistribution/stepHitType random
/gate/actor/doseDistribution/setPosition 0 0 0 cm
# fantom podzielony na macierz, nie interesuje nas pomiar wgłąb (gdyby interesował, należałoby podać też trzecią rozdzielczość)
/gate/actor/doseDistribution/setResolution 100 100 1
/gate/actor/doseDistribution/saveEveryNSeconds 60
/gate/actor/doseDistribution/enableEdep true # będziemy mierzyć zdeponowaną energię
/gate/actor/doseDistribution/enableUncertaintyEdep false # innych możliwych wielkości (niepewności, dawki, ilości zliczeń) na razie nie będziemy mierzyć
/gate/actor/doseDistribution/enableDose false
/gate/actor/doseDistribution/enableUncertaintyDose false
/gate/actor/doseDistribution/enableNumberOfHits false
/gate/actor/addActor SimulationStatisticActor stat
/gate/actor/stat/save output/stat-wynik.txt
/gate/actor/stat/saveEveryNSeconds 60
Po zdefiniowaniu geometrii, zjawisk fizycznych oraz aktora, a przed zdefiniowaniem źródeł i uruchomieniem symulacji należy dokonać inicjalizacji
#=====================================================
# INITIALISATION
#=====================================================
/gate/run/initialize
Teraz zdefiniujemy pojedyncze źródełko.
# D E F I N E T H E S O U R C E
#####
/gate/source/addSource zrodlo
/gate/source/zrodlo/gps/type Volume
/gate/source/zrodlo/gps/shape Cylinder
/gate/source/zrodlo/gps/radius 2. cm
/gate/source/zrodlo/gps/halfz 14.5 cm
/gate/source/zrodlo/gps/centre 0. 0. 30. cm
/gate/source/zrodlo/gps/particle gamma
/gate/source/zrodlo/gps/energy 140. keV # taka jest energia promieniowania Tc 99
/gate/source/zrodlo/setActivity 30000. Bq
/gate/source/zrodlo/gps/angtype iso
W niektórych wypadkach może być wygodnie podłączyć do zdefiniowanego wcześniej obiektu, na przykład tak (jeśli wcześniej zdefiniowaliśmy obiekt kolorowa_pastylka):
/gate/source/zrodlo/attachTo kolorowa_pastylka
Jeśli źródełka nie dołączyliśmy do żadnego zewnętrznego obiektu, i tak można je pokazać:
/gate/source/zrodlo/visualize
Zanim uruchomimy symulację, powinniśmy zdefiniować generator liczb pseudolosowych oraz warunek, po którym zakończy się symulacja.
# R A N D O M
# Definiujemy rodzaj generatora liczb pseudolosowych i ziarno początkowe dla jego algorytmu
# Do wyboru: JamesRandom, Ranlux64 lub MersenneTwister
/gate/random/setEngineName Ranlux64
/gate/random/setEngineSeed auto
# E X P E R I M E N T
# Definiujemy czas trwania eksperymentu
/gate/application/setTimeSlice 60. s
/gate/application/setTimeStart 0. s
/gate/application/setTimeStop 60. s
# Choć można też zamiast liczyć czasu, liczyć wyprodukowane cząstki:
# /gate/application/setTotalNumberOfPrimaries 350000
Jeśli chcemy ograniczyć ilość informacji wyrzucanych na terminal, możemy zmniejszyć do zera gadatliwość programu:
# V E R B O S I T Y czyli gadatliwość
/control/verbose 0
/run/verbose 0
/event/verbose 0
/tracking/verbose 0
/gate/random/verbose 0
/gate/application/noGlobalOutput
Dopiero w tym momencie możemy rozpocząć obliczenia:
# L E T' S R U N T H E S I M U L A T I O N !
/gate/application/startDAQ
Wynikiem symulacji będzie plik_wynikowy.img oraz plik_wynikowy.hdr. Pliki te znajdują się w jednym katalogu, pierwszy zawiera dane, drugi zawiera nagłówek, czyli opis, jak te dane powinny być interpretowane. Możemy wczytać je, obejrzeć i analizować za pomocą programu ImageJ.
Zadanie 3
Powtórzyć eksperyment Podgorsak et al, Design of X-ray targets for high energy linear accelerators in radiotherapy. Przyda się jeszcze parę praktycznych uwag:
# Filtr, czyli stożek = cone
/gate/world/daughters/name filtr
/gate/world/daughters/insert cone
/gate/filtr/geometry/setRmax1 12. cm
#/gate/filtr/geometry/setRmin1 0. cm
#/gate/filtr/geometry/setRmax2 0. cm
#/gate/filtr/geometry/setRmin2 0. cm
/gate/filtr/geometry/setHeight 3.4 cm #25.0 cm
/gate/filtr/placement/setTranslation 0. 0. 40. cm
/gate/filtr/placement/setRotationAxis 0 1 0
/gate/filtr/placement/setRotationAngle 180 deg
/gate/filtr/setMaterial Aluminium
/gate/filtr/vis/setColor yellow
#/gate/linac/vis/forceWireframe #forceSolid
#UWAGA! dodać do skryptu po inicjalizacji
#=====================================================
# ELECTRON BEAM
#=====================================================
/gate/source/addSource mybeam gps #general particle source
/gate/source/mybeam/gps/particle e-
/gate/source/mybeam/gps/pos/type Beam
/gate/source/mybeam/gps/pos/rot1 0 1 0 #wektory wyznaczające płaszczyznę przekroju wiązki
/gate/source/mybeam/gps/pos/rot2 1 0 0
/gate/source/mybeam/gps/pos/shape Circle
/gate/source/mybeam/gps/pos/centre 0 0 -3 cm #położenie
/gate/source/mybeam/gps/pos/sigma_x 1 mm
/gate/source/mybeam/gps/pos/sigma_y 1 mm
/gate/source/mybeam/gps/ene/mono 25 MeV
/gate/source/mybeam/gps/ene/type Gauss #typ rozmycia w funkcji energii
/gate/source/mybeam/gps/ene/sigma 0.5 MeV
/gate/source/mybeam/gps/direction 0 0 1
Pod adresem [6] można znaleźć różne inne typy actora, np. Energy spectrum. Daje on możliwość zapisania widma energii w formacie .root. Aby je obejrzeć należy w terminalu otworzyć program o nazwie root i napisać magiczną komendę TBrowser t (która nie jest tak bardzo magiczna, tworzy po prostu nowy obiekt o nazwie t typu TBrowser; ale roota nie będziemy się szczegółowo uczyć; aby go opuścić trzeba napisać .quit) i otworzyć plik, który nas interesuje.
Zadanie 4
Poniższy skrypt w programie root powinien narysować histogram pozycji z singli, których pozycja z była dodatnia. Analogicznie można by narysować histogram energii singli, których pozycja z była dodatnia. Itd. Można to podejście zastosować do policzenia ilości singli (lub hitów) spełniających jakiś warunek.
{
TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("kalibracja2.root");
if (!f) {
f = new TFile("kalibracja2.root");
}
//TFile *f = TFile::Open("benchSPECT.root");
TTree *Singles = (TTree*)gDirectory->Get("Singles");
//Float_t energy;
Float_t gPZ, sPZ;
//Singles->SetBranchAddress("energy",&energy);
Singles->SetBranchAddress("globalPosZ",&gPZ);
Singles->SetBranchAddress("sourcePosZ",&sPZ);
//TH1F *ener = new TH1F("ener","energia",100,0,0.2);
TH1F *posZ = new TH1F("pozycjaZ", "pozycja Z", 100,0,150);
Int_t nentries = Singles->GetEntries();
Int_t nbytes = 0;
for (Int_t i=0; i<nentries;i++) {
nbytes += Singles->GetEntry(i);
if (sPZ>0)
posZ.Fill(gPZ);
// ener.Fill(energy);
}
// **************************************** Plots **********************************************
TCanvas cont("contours","kontury",100,100,800,600);
cont.SetFillColor(0);
posZ->SetFillColor(2);
posZ->SetFillStyle(3023);
//posZ->GetXaxis()->SetTitle("");
posZ.Draw();
}