Warsztaty z Metod Terapeutycznych: Różnice pomiędzy wersjami

Z Brain-wiki
 
(Nie pokazano 8 wersji utworzonych przez 2 użytkowników)
Linia 208: Linia 208:
  
  
−
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.
+
Wynikiem symulacji będzie <tt>plik_wynikowy.img</tt> oraz <tt>plik_wynikowy.hd</tt>r. 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:
 +
 
 +
<source lang='c'>
 +
# 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
 +
 
 +
</source>
 +
 
 +
<source lang='c'>
 +
#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
 +
 
 +
</source>
 +
 
 +
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

  1. Wejdź do katalogu "Warsztaty z metod terapeutycznych" na pulpicie. Rozpakuj plik "vGATE72.vdi.7z" (jeśli nie jest już rozpakowany).
  2. 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".
  3. 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"
  4. 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/
  5. Przydatne linki
    • Strona Virtual Gate [1]
    • Wiki na temat vgate 7.2 [2]
    • Fizyka w vgate 7.2 [3]

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();

}
Źródło: „http://brain.fuw.edu.pl/edu/index.php?title=Warsztaty_z_Metod_Terapeutycznych&oldid=6279”