Luận văn Thạc sĩ Vật lý: Nghiên cứu sử dụng hiệu quả chương trình MCNP5 trong bài toán mô phỏng phổ gamma

BỘ GIÁO DỤC VÀ ĐÀO TẠO

TRƯỜNG ĐẠI HỌC SƯ PHẠM THÀNH PHỐ HỒ CHÍ MINH

TRẦN BÌNH TRANG

LUẬN VĂN THẠC SĨ VẬT LÝ

Thành phố Hồ Chí Minh – 2011

BỘ GIÁO DỤC VÀ ĐÀO TẠO

TRƯỜNG ĐẠI HỌC SƯ PHẠM THÀNH PHỐ HỒ CHÍ MINH

TRẦN BÌNH TRANG

Chuyên ngành: Vật lý nguyên tử, hạt nhân và năng lượng cao Mã số: 60. 44. 05

LUẬN VĂN THẠC SĨ VẬT LÝ

Người hướng dẫn khoa học: TS. VÕ XUÂN ÂN

Thành phố Hồ Chí Minh - 2011

LỜI CẢM ƠN

Trong quá trình thực hiện và hoàn thành luận văn này, tôi đã nhận được sự quan tâm và giúp

đỡ rất lớn từ Thầy cô và gia đình. Tôi xin được bày tỏ lòng biết ơn chân thành của mình đến:

Thầy TS. Võ Xuân Ân, người hướng dẫn khoa học, người đã mang đến cho tôi những kiến

thức phong phú, định hướng phương pháp nghiên cứu khoa học, truyền đạt tinh thần học hỏi, tìm tòi

và tận tình chỉ dẫn tôi, giúp tôi vượt qua những trở ngại, vướng mắc trong suốt quá trình thực hiện

luận văn.

Thầy TS. Nguyễn Văn Hoa, Thầy PGS. TS. Lê Văn Hoàng, hai người Thầy đã giúp đỡ tôi

vượt qua giai đoạn khó khăn, đóng góp ý kiến và động viên tôi từ những ngày đầu thực hiện luận

văn.

Thầy TS. Nguyễn Mạnh Hùng, Trưởng khoa Vật lý, trường Đại học Sư Phạm TPHCM đã tạo

điều kiện cho tôi được thực hiện luận văn này.

Cuối cùng, xin cảm ơn gia đình tôi, bạn bè tôi đã luôn bên cạnh động viên và giúp đỡ tôi

trong những lúc khó khăn nhất.

MỤC LỤC

BẢNG CÁC CHỮ VIẾT TẮT

Chữ viết tắt Tiếng Việt Tiếng Anh

ACTL Thư viện số liệu ACTL ACTivation Library

ENDF Thư viện số liệu ENDF Evaluated Nuclear Data File

ENDL Thư viện số liệu ENDL Evaluated Nuclear Data Library

Ge Germanium Germanium

Detector germanium siêu tinh khiết High Purity Germanium HPGe

KHTN Khoa học tự nhiên -

MCA Multi Channel Analyzer

MCNP Monte Carlo N – Particle

NAS Khối phân tích biên độ đa kênh Chương trình mô phỏng Monte Carlo MCNP Hãng cung cấp nguồn phóng xạ North American Scientific, Inc.

NEA Nuclear Energy Agency

NJOY - Cơ quan năng lượng hạt nhân Mã định dạng các thư viện số liệu hạt nhân trong MCNP

NSS Hãng cung cấp nguồn phóng xạ Nuclear Services & Supplies Rost GmbH

TP HCM Thành phố Hồ Chí Minh -

MỞ ĐẦU

Vật lý hạt nhân ngày nay đang có những bước phát triển rất mạnh mẽ, đặc biệt trong lĩnh vực

khoa học hạt nhân ứng dụng. Từ khi phát hiện ra hiện tượng phóng xạ, việc nghiên cứu các hiện

tượng vật lý hạt nhân dựa trên đo đạc phổ phóng xạ ngày càng trở nên phổ biến. Trong đó, lĩnh vực

đo phổ gamma được tập trung nghiên cứu và đem lại nhiều kết quả thực tiễn quan trọng. Hiện nay,

công nghệ đo phổ gamma được phát triển ở mức độ cao và được sử dụng phổ biến trong các phòng

nghiên cứu. Tại Việt Nam, Viện Khoa học và Kỹ thuật Hạt nhân Hà Nội, Viện Nghiên cứu Hạt nhân

Đà Lạt, Trung tâm Hạt nhân Thành phố Hồ Chí Minh (TPHCM) và một số phòng thí nghiệm vật lý

hạt nhân thuộc các trường đại học đã được trang bị các hệ phổ kế gamma để phục vụ nghiên cứu và

đo đạc, khảo sát các mẫu môi trường.

Có hai cách để tiến hành khảo sát phổ gamma, đo đạc trực tiếp hoặc sử dụng chương trình

mô phỏng. Một trong những chương trình mô phỏng đang được sử dụng rộng rãi hiện nay để giải

quyết các vấn đề trong vật lý hạt nhân là chương trình MCNP5. Đây là một chương trình mô phỏng

có độ tin cậy cao vì đã được kiểm chứng và sử dụng trong nhiều năm qua và ở nhiều phòng thí

nghiệm trong nước cũng như trên thế giới.

Đã có nhiều nghiên cứu trong nước và trên thế giới sử dụng chương trình MCNP5 để xây

dựng đường cong hiệu suất, phân bố suất liều, nghiên cứu lò phản ứng … nhưng tiến hành riêng lẻ

(chạy đơn). Một số công trình nghiên cứu trong nước gần đây như “tính toán tối ưu hộp chứa mẫu

phóng xạ” [8], “xây dựng ma trận đáp ứng nhờ phổ gamma đơn năng” [5] và nhiều công trình

nghiên cứu khác đã sử dụng chương trình MCNP5 mô phỏng phổ gamma để giải quyết bài toán.

Tuy nhiên, khó khăn chính mà các tác giả này gặp phải đó là trong quá trình mô phỏng phải tính

MCNP5 liên tục và nhiều lần, trong khi vẫn chưa có một công cụ giúp kết nối giữa MCNP5 và

chương trình tính toán. Do đó công việc buộc phải thực hiện rời rạc và vì vậy tốn kém nhiều thời

gian để thiết lập input và triển khai chạy chương trình MCNP5.

Việc nghiên cứu để kết nối giữa MCNP5 với chương trình tính toán sẽ giúp gia tăng đáng kể

hiệu quả của chương trình MCNP5 trong bài toán mô phỏng phổ gamma, nhờ đó rút ngắn thời gian

và giảm bớt khối lượng công việc cho các nhà nghiên cứu. Không những vậy, tự động hóa quá trình

tính toán giúp giảm thiểu sai sót, đảm bảo tính chính xác của kết quả.

Việc xây dựng chương trình tính toán kết nối với MCNP5 và tích hợp các modun tính toán

xử lý phổ gamma sẽ hỗ trợ người dùng MCNP5 trong việc giải quyết các bài toán vật lý hạt nhân,

giảm thiểu được thời gian và công sức thực hiện các quá trình tính toán lặp lại. Nhờ đó chương trình

MCNP5 sẽ được khai thác hiệu quả hơn, đặc biệt đối với những bài toán phải mô phỏng một số

lượng lớn phổ gamma.

Từ những phân tích trên tôi đã chọn đề tài: “Nghiên cứu sử dụng hiệu quả chương trình

MCNP5 trong bài toán mô phỏng phổ gamma” làm luận án tốt nghiệp.

Mục tiêu của luận án là: (1) Dựa trên phương pháp mô phỏng Monte Carlo để xây dựng các

modun tính toán trong bài toán mô phỏng phổ gamma; (2) Thiết lập chương trình tính toán bằng

ngôn ngữ lập trình Fortran; (3) Kết nối chương trình MCNP5 và chương trình tính toán, đảm bảo

quá trình tính toán được liền mạch, không rời rạc, đơn lẻ; (4) Đánh giá về kết quả thu được và đưa

ra các kiến nghị liên quan.

Đối tượng nghiên cứu của luận án này là chương trình MCNP5, đây là chương trình mô

phỏng do Phòng thí nghiệm quốc gia Los Alamos, Hoa Kỳ xây dựng. MCNP5 là một công cụ tính

toán mạnh, có thể mô phỏng các quá trình vận chuyển của neutron, photon, electron trong môi

trường vật chất. Phần mềm này được xây dựng dựa trên phương pháp mô phỏng Monte Carlo.

Phương pháp nghiên cứu của luận án là sử dụng ngôn ngữ lập trình Fortran để xây dựng

chương trình tính toán kết nối với MCNP5 trong bài toán mô phỏng phổ gamma.

Với nội dung đó, luận văn được trình bày thành bốn phần như sau:

+ Chương 1: TỔNG QUAN, giới thiệu một cách khái quát các vấn đề về tương tác của

photon với vật chất, về phương pháp mô phỏng Monte Carlo với chương trình MCNP5, về ngôn

ngữ lập trình Fortran.

+ Chương 2: MÔ HÌNH HÓA HỆ PHỔ KẾ GAMMA BẰNG CHƯƠNG TRÌNH MCNP5,

mô tả hệ phổ kế gamma, trình bày cách thiết lập input và các bước thực hiện một bài toán mô phỏng

bằng chương trình MCNP5.

+ Chương 3: XÂY DỰNG CHƯƠNG TRÌNH TÍNH TOÁN, xây dựng các modun tính toán

để sử dụng hiệu quả chương trình MCNP5 trong bài toán mô phỏng phổ gamma bằng ngôn ngữ lập

trình Fortran.

+ Chương 4: KẾT LUẬN VÀ KIẾN NGHỊ, tổng kết và đánh giá các kết quả đạt được, đưa ra

kiến nghị về những hướng nghiên cứu tiếp theo.

CHƯƠNG 1: TỔNG QUAN

1.1. TƯƠNG TÁC GIỮA PHOTON VỚI MÔI TRƯỜNG VẬT CHẤT – PHỔ

GAMMA

1.1.1. Tương tác của bức xạ gamma với vật chất

Bức xạ gamma tương tác với môi trường vật chất thông qua các quá trình hấp thụ và tán xạ.

Trong quá trình hấp thụ, tia gamma sẽ truyền toàn bộ năng lượng cho các hạt trong vật chất, sau đó

tia gamma biến mất. Còn trong quá trình tán xạ, tia gamma chỉ truyền cho các hạt vật chất một phần

năng lượng và nó bị tán xạ dưới một góc nào đó (phương chuyển động ban đầu bị thay đổi). Quá

trình tương tác giữa tia gamma và vật chất được gọi là sự ion hóa gián tiếp vì các sản phẩm tạo ra

sau va chạm (các hạt vi mô tích điện hay các photon thứ cấp) sẽ tác dụng tiếp với các hạt trong môi

trường vật chất và tạo ra phần lớn các ion.

Các tia gamma có thể tương tác với vật chất theo nhiều cơ chế khác nhau, tuy nhiên, trong

ghi đo phóng xạ, người ta dựa vào ba quá trình đóng vai trò quan trọng nhất: hấp thụ quang điện, tán

xạ Compton và hiệu ứng tạo cặp [15]. Một số hiệu ứng khác như tán xạ Thomson, phản ứng quang

hạt nhân,… có xác suất thấp nên có thể bỏ qua.

1.1.1.1. Hấp thụ quang điện

Hấp thụ quang điện là hiện tượng photon tới bị electron trong nguyên tử hấp thụ hoàn toàn

năng lượng, sau đó photon biến mất, electron bứt ra khỏi lớp vỏ nguyên tử. Các electron này được

gọi là electron quang điện.

lkε , tùy thuộc quỹ đạo

Mỗi electron quỹ đạo ứng với một giá trị năng lượng liên kết xác định

chuyển động (K, L, M, N,…) và số nguyên tử Z của hạt nhân. Năng lượng liên kết của nguyên tử

lớn hơn đối với các electron ở vỏ nguyên tử sâu hơn và với nguyên tử có số Z lớn hơn. Hiệu ứng

quang điện mạnh nhất xảy ra đối với lượng tử gamma có năng lượng có thể so sánh được với năng

lượng liên kết của nguyên tử.

lkε của electron để có thể phá vỡ liên kết của

Năng lượng của photon tới E phải lớn hơn

electron với hạt nhân. Phần năng lượng dư thừa chính là động năng cho quang electron:

lkε + ERe R (1.1)

E = hν =

Với động năng EReR đó, quang electron có khả năng ion hóa các nguyên tử và phân tử khác.

Phần động năng EReR của quang electron lớn hơn rất nhiều so với phần năng lượng để bứt electron ra

lkε .

khỏi quỹ đạo

Về phía nguyên tử vật chất, khi một electron bị bật ra khỏi quỹ đạo, electron khác ở vành

lkε giữa hai quỹ đạo sẽ được

ngoài có thể đến thế chỗ. Năng lượng dư thừa do sự chênh lệch của

lkε phụ thuộc vào các quỹ đạo và vào số nguyên

phát ra dưới dạng một photon. Giá trị năng lượng

tử Z nên photon thứ cấp này có giá trị xác định và được gọi là bức xạ đặc trưng.

lkε = Kε thì hiệu ứng quang điện không xảy ra với các electron lớp K mà chỉ xảy ra

Nếu E <

lkε =

Lε thì hiệu ứng quang điện không xảy ra với các

với các electron ở lớp L, M, N,… nếu E <

Lε > Mε ).

electron lớp K, L mà chỉ xảy ra với các electron ở lớp M, N, O,… (Vì Kε >

Hiệu ứng quang điện không xảy ra với electron tự do vì không bảo đảm quy luật bảo toàn

năng lượng và động lượng.

Như vậy muốn có hiệu ứng quang điện thì:

Electron phải liên kết trong nguyên tử •

Năng lượng tia gamma phải lớn hơn năng lượng liên kết của electron nhưng không •

được lớn quá vì khi đó electron có thể coi gần như tự do.

Nhận xét này được thể hiện trên hình mô tả sự phụ thuộc tiết diện hiệu ứng quang điện vào

năng lượng gamma

1/E7/2

0 E

Hình 1.1: Tiết diện hiệu ứng quang điện phụ thuộc năng lượng gamma E

Ở miền năng lượng gamma lớn thì tiết diện rất bé vì khi đó gamma coi electron là liên kết rất

1 E

ε K E

yếu. Khi giảm năng lượng gamma, tức là tăng tỉ số , tiết diện tăng theo quy luật . Khi E tiến

1 7/2 E

dần đến Kε , tiết diện tăng theo hàm và tăng cho đến khi E = Kε . Khi năng lượng gamma vừa

giảm xuống dưới giá trị Kε thì hiệu ứng quang điện không thể xảy ra với electron lớp K nữa nên tiết

diện giảm đột ngột. Tiếp tục giảm năng lượng gamma, tiết diện tăng trở lại do hiệu ứng quang điện

Lε rồi lại giảm đột ngột khi E giảm xuống thấp hơn

Lε . Sau đó hiệu ứng quang điện xảy ra đối với electron lớp M,…

đối với electron lớp L. Nó đạt giá trị lớn tại E =

5 thuộc vào Z, theo quy luật ZP

P.

Do năng lượng liên kết thay đổi theo số nguyên tử Z nên tiết diện tương tác quang điện phụ

Sự đóng góp của hiệu ứng quang điện đối với các lớp L, M,… bé so với electron lớp K. Các

nghiên cứu cho thấy hiệu ứng quang điện xảy ra chủ yếu với electron lớp K và với tiết diện rất lớn

đối với các nguyên tử nặng, chẳng hạn chì, ngay cả ở vùng năng lượng cao. Còn đối với các nguyên

tử nhẹ, chẳng hạn cơ thể sinh học, hiệu ứng quang điện chỉ xuất hiện ở vùng năng lượng thấp.

1.1.1.2. Hiệu ứng Compton

Khi tăng năng lượng gamma đến giá trị lớn hơn nhiều so với năng lượng liên kết của electron

K trong nguyên tử thì vai trò của hiệu ứng quang điện không còn đáng kể và bắt đầu hiệu ứng

Compton. Khi đó có thể bỏ qua năng lượng liên kết của electron so với năng lượng gamma và tán xạ

gamma lên electron có thể coi như tán xạ với electron tự do, gọi là tán xạ Compton.

Hiệu ứng Compton là sự tán xạ đàn hồi của gamma vào các electron chủ yếu ở quỹ đạo ngoài

cùng của nguyên tử. Sau quá trình tán xạ, lượng tử gamma thay đổi phương bay và bị mất một phần

năng lượng còn electron được giải phóng ra khỏi nguyên tử.

Hình 1.2: Hiệu ứng Compton

Trên cơ sở tính toán động học quá trình tán xạ đàn hồi của hạt gamma chuyển động với năng

lượng E lên electron đứng yên. Gọi năng lượng gamma sau tán xạ là E’, năng lượng electron sau tán

xạ EReR và θ là góc bay của gamma sau tán xạ, ta có

=

E

'

E

1 −

+

α

1

θ (1 cos )

=

E

(1.2)

eE

α +

− θ (1 cos ) − θ α (1 cos )

1

-31

2

α=

(1.3)

8 P m/s, P kg, c = 3.10P

em = 9,1.10P

em c = 0,51 MeV

2

E m c e

, Trong đó:

1

φ= −

Góc bay φ của electron sau tán xạ liên hệ với góc θ như sau:

tg

ctg

θ 2

−

1

E E

'

(1.4)

’ 'λ của gamma liên hệ với các giá trị năng lượng E và EP P của nó như sau:

λ=

λ = '

Các bước sóng λ và

hc E

hc E '

’ Theo công thức (1.2) thì EP P < E, nghĩa là năng lượng gamma giảm sau tán xạ Compton và

; (1.5)

bước sóng của nó tăng. Gia số tăng bước sóng phụ thuộc vào góc tán xạ θ của gamma theo biểu

∆ =

θ

− = '

2

2 sin (

/ 2)

thức:

λ λ λ λ c

−

12

=

=

m

2, 42.10

(1.6)

λ c

h m c e

λ∆ chỉ phụ thuộc vào góc θ nên không phụ thuộc vào vật liệu của môi trường. Từ công thức (1.6),

θ π=

/ 2

'λ tăng khi tăng góc tán xạ và

0λ∆ = khi

0ϕ= ;

là bước sóng Compton, được xác nhận bởi thực nghiệm. Do Trong đó:

cλ λ∆ =

0λ∆ = khi θ π= ; Tuy nhiên với một góc θ cho trước thì λ∆ không phụ thuộc vàoλ. Như vậy hiệu

∆ << vì khi đó

'λ λ= , chẳng hạn với ánh sáng

khi ; có thể thấy rằng bước sóng

ứng Compton không đóng vai trò đáng kể khi λ λ

nhìn thấy hoặc ngay cả với tia X năng lượng thấp. Hiệu ứng Compton chỉ đóng góp lớn đối với tia

gamma sóng ngắn, hay năng lượng cao, sao cho λ λ∆ ≈ .

Góc bay θ của gamma tán xạ có thể thay đổi từ 0o đến 180o. Khi tán xạ Compton, năng

lượng tia gamma giảm và phần năng lượng giảm đó truyền cho electron giật lùi. Như vậy năng

θ=

180

o, tức là khi tán xạ giật lùi (hay còn gọi là tán xạ

lượng electron giật lùi càng lớn khi gamma tán xạ với góc θ càng lớn. Gamma truyền năng lượng

lớn nhất cho electron khi tán xạ ở góc

ngược).

2

2

θ

=

Tiết diện vi phân của tán xạ Compton có dạng:

2 r e

2 α 2

σ d Ω d

) θ

+

(1.7)

( + 1 cos

+ 1 cos ( α − 1 cos

) θ

) ( θ − 1 cos ) ( − + θ α  1 1 cos 

 

    

  + 1 2    

 2 1 

2

=

=

α ;

r e

2

2

e m c e

E m c e

Trong đó

Tiết diện tán xạ Compton toàn phần nhận được bằng cách lấy tích phân biểu thức (1.7) theo

+

=

+

+

+

−

σ

tất cả các góc tán xạ:

( ln 1 2

) α

( ln 1 2

) α

Compton

2 π 2 er

2

) α 1 − α α

( 2 1 + 1 2

1 α 2

+ 1 3 ( + 1 2

α ) α

  

  

 + α  1  2 α  

    

(1.8)

1.1.1.3. Hiệu ứng tạo cặp

P). Đó là hiệu ứng

+ - P) và positron (eP với trường hạt nhân đó và biến chuyển thành một cặp electron (eP

Những photon có năng lượng E Rγ R≥ 1,022 MeV khi đến gần hạt nhân nguyên tử sẽ tương tác

P = 0,511 MeV).

2 khối lượng tĩnh mReR của hai hạt vi mô đó (E = mRe RcP

tạo cặp electron-positron. Năng lượng tối thiểu dùng cho hiệu ứng này là 1,022 MeV tương ứng với

Phần năng lượng còn lại của photon tới trở thành động năng của hai hạt vi mô mới xuất hiện.

Như vậy:

2 ERγ R = 2mReRcP

+ - PRd R + EP P + EP

PRd R

(1.9)

Các hạt thứ cấp này có động năng nên sẽ tương tác với vật chất và cũng gây ra quá trình ion

hóa thứ cấp. Electron sẽ mất dần động năng rồi chuyển về dạng chuyển động nhiệt hoặc gắn với một

ion dương nào đó. Positron mang điện tích dương nên sẽ dễ dàng kết hợp với các electron khác

trong vật chất, điện tích của chúng bị trung hòa, chúng hủy lẫn nhau, gọi là hiện tượng hủy electron-

positron. Khi hủy electron-positron, hai lượng tử gamma sinh ra bay ngược chiều nhau, mỗi lượng

tử có năng lượng 0,511 MeV, tức là năng lượng tổng cộng của chúng bằng tổng khối lượng hai hạt

electron và positron 1,022 MeV.

Sự biến đổi năng lượng thành khối lượng như trên phải xảy ra gần một hạt nào đó để hạt này

chuyển động giật lùi giúp tổng động lượng được bảo toàn. Quá trình tạo cặp cũng có thể xảy ra gần

electron nhưng xác suất bé so với quá trình tạo cặp gần hạt nhân khoảng 1000 lần.

Tiết diện hiệu ứng tạo cặp trong trường hạt nhân có dạng phức tạp. Sau đây là các biểu thức

2

=

−

σ

ln

tiết diện trong vài miền năng lượng của tia gamma:

pair

2 r e

Z 137

28 9

218 27

E 2 2 m c e

  

  

2

2

1/ 3

<<

<<

E

137

m c Z −

(1.10)

m c e

e

2

−

1/3

=

−

σ

Z

đối với khi không tính đến hiệu ứng màng che chắn.

pair

2 r e

Và: (1.11)

( ln 183

)

Z 137

28 9

2 27

  

  

2

1/ 3

>>

E

m c Z −

137 e

2

1/3

137

30

đối với , khi tính đến hiệu ứng màng che chắn toàn phần.

em c Z − = MeV đối với nhôm và 15MeV đối với chì.

2

2

<

<

E

5

50

2 tiết diện tạo cặp tỉ lệ với ZP

P và lnE.

Trong đó

m c e

m c e

2 Tiết diện tạo cặp electron-positron gần tỉ lệ với ZP

P nên có giá trị lớn đối với chất hấp thụ có

Trong miền năng lượng

số nguyên tử lớn.

1.1.1.4. Tổng hợp các hiệu ứng khi gamma tương tác với vật chất

Như đã trình bày ở trên, khi gamma tương tác với vật chất có ba hiệu ứng chính xảy ra, đó là

hiệu ứng quang điện, hiệu ứng Compton và hiệu ứng tạo cặp electron-positron. Tiết diện vi phân

= σ σ

+

σ

+

σ

tương tác tổng cộng của các quá trình này bằng:

photo

Compton

pair

5

(1.12)

photoσ tỷ lệ với

Z 7 / 2 E

Trong đó tiết diện quá trình quang điện

Comptonσ

Z E

2 lnZ

tỷ lệ với Tiết diện quá trình tán xạ Compton

E

pairσ tỷ lệ với

Tiết diện quá trình tạo cặp

Hình 1.3: Sự phụ thuộc năng lượng của các quá trình tương tác gamma khác nhau trong NaI

(Theo The Atomic Nuclear, R.D.Evans, 1955) [17]

Hình 1.3 là đồ thị biểu diễn sự phụ thuộc của tiết diện hấp thụ quang điện trong NaI theo

năng lượng. Ở vùng năng lượng thấp (tương ứng với năng lượng liên kết của các electron ở các lớp

khác nhau) có những mép hấp thụ xuất hiện. Những mép với năng lượng cao nhất sẽ ứng với các

electron lớp K. Ở phía trên ngay sát mép này, năng lượng photon chỉ đủ để chịu một hấp thụ quang

điện trong đó một electron lớp K bị bứt ra khỏi nguyên tử. Ở phía dưới ngay sát mép này, không có

đủ năng lượng để quá trình này xảy ra nên xác suất tương tác giảm nhanh đột ngột. Tương tự, các

mép hấp thụ ở năng lượng thấp hơn ứng với electron các lớp L, M,…trong nguyên tử.

Từ sự phụ thuộc các tiết diện vào năng lượng E của gamma và điện tích Z của vật chất như

trên, có thể thấy rằng trong miền năng lượng bé hơn ER1 R cơ chế cơ bản trong tương tác gamma với

vật chất là hiệu ứng quang điện, trong miền năng lượng trung gian ER1 R < E < E R2 R là hiệu ứng

Compton và trong miền năng lượng cao E > ER2 R là quá trình tạo cặp electron-positron. Các giá trị

năng lượng phân giới E R1 R và E R2 R phụ thuộc vào vật chất. Đối với nhôm thì ER1 R = 50 keV, ER2 R = 15

MeV. Còn đối với chì ER1 R = 500 keV và ER2 R = 5 MeV.

1.1.2. Nguyên tắc của sự hình thành phổ gamma

Vấn đề phân tích hạt nhân ngày càng trở thành nhu cầu bức thiết trong các nghiên cứu thực

nghiệm. Việc chế tạo thành công các loại detector bán dẫn germanium siêu tinh khiết (high purity

germanium – HPGe) với độ phân giải và hiệu suất đếm cao vào những năm 1980 là một bước ngoặt

trong lịch sử phát triển các thiết bị ghi nhận bức xạ tia X và tia gamma vì nó đã cải thiện đáng kể độ

chính xác của các phép phân tích bằng phương pháp hạt nhân. Sử dụng những detector ghi bức xạ

gamma trong phân tích phóng xạ có ưu điểm quan trọng là kỹ thuật phân tích không phá mẫu và khả

năng phân giải cao. Các thiết bị ghi nhận tia gamma và dựa vào những dữ liệu đó để xây dựng phổ

gamma. Phổ gamma chứa đựng những thông tin về những hiệu ứng xảy ra trong vật chất khi chiếu

bức xạ gamma vào, từ đó có thể phân tích, xác định năng lượng bức xạ gamma phát ra và hàm

lượng của các hạt nhân phóng xạ trong mẫu.

Khi đi qua môi trường của detector với cấu hình cụ thể, do bức xạ gamma không mang điện

tích nên không gây hiệu ứng ion hóa hoặc kích thích trực tiếp vào detector. Vì vậy, việc ghi nhận

chúng được thực hiện thông qua các tương tác mà trong đó một phần hoặc toàn bộ năng lượng của

chúng được truyền cho electron. Chính các electron này gây ion hóa tạo ra các xung điện ở lối ra

của detector. Như vậy detector đã biến đổi năng lượng tia gamma thành năng lượng của các electron

nhanh đồng thời ghi nhận chuyển đổi electron nhanh thành những tín hiệu điện.

Detector hoạt động theo kiểu xung, mỗi xung riêng rẽ mang thông tin quan trọng liên quan

đến điện tích được tạo ra bởi tương tác của bức xạ trong detector. Những xung này được tập hợp và

lưu trữ cho sự thể hiện phân bố biên độ xung của detector ở đầu ra. Một cách tổng quát phổ gamma

đo được là kết quả của sự tương tác của hệ detector lên phổ gamma tới, làm phân bố lại dạng của

phổ gamma tới, bao gồm đỉnh toàn phần do hiệu ứng quang điện và nền liên tục từ hiệu ứng tán xạ

Compton nhiều lần trong môi trường detector và các vật liệu xung quanh. Có hai cách thông thường

để trình bày thông tin về phân bố biên độ xung là phổ vi phân và phổ tích phân.

Quá trình này có thể cụ thể hóa bằng các bước sau:

1) Tia gamma có năng lượng hν đi vào detector Ge.

2) Các electron thứ cấp với tổng động năng E (eV) sinh ra trong vùng nhạy trong tinh thể Ge

bởi sự tương tác photon với Ge.

3) Các electron thứ cấp sinh ra một lượng lớn cặp electron – lỗ trống bởi sự ion hoá trong

tinh thể Ge.

4) Một tín hiệu xung có biên độ V (Volt) được sinh ra ở đầu vào tiền khuếch đại với điện

dung C (Farad).

5) Tín hiệu xung có biên độ V được khuếch đại và biến đổi thành số chỉ vị trí kênh (ch) bởi

khối xử lý biên độ đa kênh (MCA), một số đếm được cộng vào số đếm tổng tại vị trí kênh ch.

6) Nhiều tia gamma được phân tích bởi quá trình như trên, và sự phân bố độ cao của xung

được hình thành (phổ gamma).

Tia gamma

Electron thứ cấp

Cặp electron – lỗ trống

Tín hiệu xung

Phân tích độ cao xung

Phổ gamma

Hình 1.4: Sơ đồ sự hình thành phổ gamma

Như đã nói ở trên, phổ gamma đo được là kết quả của sự tương tác của hệ detector lên phổ

gamma tới, làm phân bố lại dạng của phổ gamma tới, sự phân bố này do hiệu ứng quang điện và

hiệu ứng tán xạ Compton nhiều lần trong môi trường detector và các vật liệu xung quanh. Kết quả

của hấp thụ quang điện là giải phóng các electron quang điện (mang hầu hết năng lượng của

gamma) cùng với một hoặc một số electron năng lượng thấp hơn ứng với sự hấp thụ năng lượng liên

kết của electron quang điện. Nếu không có sự thất thoát ra khỏi detector thì tổng động năng của các

electron được tạo ra phải bằng với năng lượng ban đầu của photon. Vì thế hấp thụ quang điện là một

quá trình lý tưởng cho việc đo đạc năng lượng của gamma. Với chùm gamma đơn năng và những

điều kiện lý tưởng, tổng động năng của các electron bằng với năng lượng gamma tới và phân bố vi

phân của động năng electron sau mỗi chuỗi các sự kiện hấp thụ quang điện sẽ có dạng một hàm

delta đơn giản như hình 1.5. Một đỉnh đơn xuất hiện tại năng lượng ứng với năng lượng của gamma

tới.

Hình 1.5: Đỉnh năng lượng toàn phần trong phổ độ cao xung vi phân

Trong tán xạ Compton, thông thường gamma có thể bị tán xạ ở bất kỳ góc nào khi tương tác

xảy ra trong detector. Do đó electron cũng có thể nhận một năng lượng bất kỳ từ 0 đến giá trị cực

đại ứng với θ π= và phân bố năng lượng electron có dạng tổng quát như hình dưới.

Hình 1.6: Nền Compton trong phổ độ cao xung vi phân

Trong quá trình tạo cặp, tổng động năng của các hạt tích điện (electron và positron) được tạo

bởi gamma tới cũng có dạng hàm delta đơn giản.

Hình 1.7: Đỉnh tạo cặp trong phổ độ cao xung vi phân

Tạo cặp là một quá trình phức tạp bởi positron là hạt không bền và chỉ đi được vài milimet.

Khi bị làm chậm trong môi trường hấp thụ đến năng lượng cỡ năng lượng nhiệt của electron,

positron sẽ hủy với một electron và một cặp photon 0,511 MeV xuất hiện. Có ba khả năng xảy ra:

hν MeV. Như vậy tia gamma mất hoàn toàn năng lượng nên có sự đóng góp vào số đếm toàn phần.

+ Cả hai photon đều bị hấp thụ. Năng lượng của tia gamma bị mất là: ( hν-1,022 + 1,022) =

+ Chỉ có một photon hủy bị hấp thụ, một photon thoát ra ngoài nên năng lượng tia gamma

mất trong vùng nhạy là: hν- 1,022 + 0,511 = hν- 0,511 MeV. Các xung này đóng góp số đếm vào

phổ biên độ xung tạo thành đỉnh thoát cặp thứ nhất (đỉnh thoát đơn, SE).

+ Khi cả hai photon hủy đều thoát khỏi tinh thể, năng lượng tia gamma mất trong vùng nhạy

của detector là: hν- 1,022 MeV. Các xung này đóng góp số đếm vào phổ biên độ xung tạo thành

đỉnh thoát cặp thứ hai (đỉnh thoát đôi, DE).

Tuy nhiên, các xác suất SE và DE này thường rất thấp. Do đó với nguồn cường độ mạnh

hoặc đo thời gian dài mới khảo sát được các đỉnh này.

Thời gian của quá trình hủy electron và positron rất ngắn, vì thế bức xạ hủy xuất hiện gần

như cùng lúc với tương tác tạo cặp ban đầu. Bức xạ hủy này gây ảnh hưởng đáng kể lên đáp ứng của

1.2. PHƯƠNG PHÁP MÔ PHỎNG MONTE CARLO VÀ CHƯƠNG TRÌNH MCNP5

1.2.1. Phương pháp mô phỏng Monte Carlo

các detector gamma.

1.2.1.1. Giới thiệu về mô phỏng

• Giới thiệu

Trên thế giới, trong nhiều lĩnh vực khoa học, để đáp ứng nhu cầu tìm hiểu và khám phá của

các nhà nghiên cứu, cùng với sự xuất hiện của máy tính thì việc mô phỏng các hiện tượng thí

nghiệm trở nên hết sức cần thiết.

Mô phỏng là việc sử dụng máy tính kết hợp với các quy luật toán học, vật lý, dựa trên phép

tính định lượng, tương đối hoá các tham số để giải các bài toán, nghiên cứu kết cấu hay quá trình,

thực hiện tính toán hay dựng lên các mô hình của các thí nghiệm.

• Các loại mô phỏng

Để thực hiện một bài toán mô phỏng dù đơn giản hay phức tạp thì đều phải mô hình hóa và

lựa chọn phương thức thích hợp để thực hiện trên máy tính.

Dựa trên nhiều tiêu chuẩn, mô phỏng có thể chia thành nhiều loại. Dưới đây là vài loại cơ

bản:

- Mô phỏng ngẫu nhiên: còn gọi là mô phỏng Monte Carlo, nó áp dụng nguyên tắc gieo các

số ngẫu nhiên để mô phỏng các hiện tượng ngẫu nhiên.

- Mô phỏng tất định: là phương pháp tính toán có thể đoán trước được. Nếu nó chạy với một

dữ liệu vào cụ thể thì các dữ liệu ra không đổi.

- Mô phỏng liên tục: bằng việc sử dụng các phương trình vi phân và giải tích số, máy tính sẽ

giải phương trình một cách tuần hoàn và sử dụng kết quả thu được để thay đổi trạng thái, số liệu

xuất ra.

- Mô phỏng rời rạc: người ta ghi lại một dãy các sự kiện đã được sắp xếp theo thời gian, khi

mô phỏng các sự kiện này sẽ tạo ra các sự kiện mới.

1.2.1.2. Phương pháp mô phỏng Monte Carlo

Thuật ngữ “phương pháp Monte Carlo” xuất hiện từ Thế chiến thứ hai khi tiến hành các mô

phỏng ngẫu nhiên trong quá trình chế tạo bom nguyên tử. Phương pháp Monte Carlo là phương

pháp giải số, áp dụng cho việc mô phỏng sự tương tác giữa các vật thể hay giữa vật thể với môi

trường nhờ lý thuyết cơ học và động lực học, dựa theo yêu cầu của hệ cần mô phỏng.

Đây là một phương pháp ngẫu nhiên, nhờ sự phát sinh các số ngẫu nhiên để tính toán.

Phương pháp này thường được dùng để mô phỏng các quá trình ngẫu nhiên của hệ, nó không liên

quan đến thời gian nên không thể sử dụng để mô phỏng các đại lượng phụ thuộc thời gian.

Các phương pháp Monte Carlo rất khác với các phương pháp vận chuyển tất định. Các

phương pháp tất định giải phương trình vận chuyển đối với trạng thái hạt trung bình. Ngược lại,

phương pháp Monte Carlo không giải phương trình tường minh mà là sự mô phỏng các hạt riêng rẽ

và ghi lại một số khía cạnh trạng thái trung bình của chúng. Trạng thái trung bình của các hạt trong

vật lý khi đó được rút ra từ trạng thái trung bình của các hạt được mô phỏng (bằng cách sử dụng

định lý giới hạn trung tâm). Phương pháp Monte Carlo cho phép biểu diễn chi tiết tất cả các khía

cạnh của các số liệu vật lý trong quá trình vận chuyển hạt.

Trong những năm gần đây, các chương trình mô phỏng vận chuyển bức xạ bằng phương

pháp Monte Carlo được sử dụng ngày càng rộng rãi. Điều này, một mặt do yêu cầu cấp bách giải

quyết nhiều bài toán quan trọng trong thực tế từ thiết kế lò phản ứng đến bảo vệ bức xạ và vật lý y

học. Mặt khác, các chi phí thực nghiệm tăng lên và các chi phí tính toán giảm xuống cũng làm cho

việc mô phỏng trở nên hấp dẫn hơn, đặc biệt khi các thí nghiệm được tiến hành trong môi trường

nguy hiểm. Hơn nữa, các kỹ thuật tính trong các chương trình này cũng nhanh hơn và các máy tính

cũng tốt hơn làm cho mô phỏng số trở nên đáng tin cậy.

Đa số các công trình nghiên cứu trên thế giới tập trung vào các vấn đề về mô phỏng đáp ứng

phổ, sử dụng mô phỏng trong việc hỗ trợ tính toán hiệu suất đối với các dạng hình học nguồn và

mẫu khác nhau, khảo sát hiệu suất theo năng lượng, theo khoảng cách, hiệu chỉnh trùng phùng tổng

đối với gamma phân rã nhiều tầng, hiệu chỉnh tự hấp thụ đối với hình học nguồn và mẫu thể tích,…

Phương pháp Monte Carlo phụ thuộc rất nhiều vào tập hợp các số ngẫu nhiên. Trên máy tính,

các số ngẫu nhiên trong các ngôn ngữ lập trình thường chỉ có giá trị từ 0 đến 1, trong khi phạm vi

của các đại lượng cần mô phỏng thì rất rộng lớn. Vì vậy cần có những kĩ thuật thích hợp để mô

phỏng bài toán.

Phương pháp này được ứng dụng rất rộng rãi trong nhiều lĩnh vực như thống kê, lượng tử,

mô phỏng cấu trúc vật liệu và cả trong kinh tế nữa. Trên thực tế, nó được sử dụng trong nghiên cứu

1.2.2. Chương trình MCNP5

sự vận chuyển bức xạ, tính toán và thiết kế lò phản ứng …

1.2.2.1. Tổng quan về MCNP

MCNP (Monte Carlo N-Particle) là chương trình phổ biến để mô phỏng tương tác giữa

neutron, gamma, electron với môi trường vật chất.

MCNP ban đầu được phát triển bởi nhóm Monte Carlo và sau này bởi nhóm Radiation

Transport (Nhóm X-5) của phòng Vật lý Lý thuyết Ứng dụng ở Phòng thí nghiệm quốc gia Los

Alamos (Hoa Kỳ). Nhóm X-5 cải tiến MCNP và cứ hai hoặc ba năm họ lại cho ra một phiên bản

MCNP mới, các phiên bản này được cung cấp tới người dùng thông qua Trung tâm Thông tin Che

chắn Bức xạ (RSICC) ở Oak Ridge, Tennessee (Hoa Kỳ) và ngân hàng dữ liệu OECD/NEA ở Pari

(Pháp). MCNP sử dụng các thư viện số liệu hạt nhân và nguyên tử từ các nguồn dữ liệu ENDF

(Evaluated Nuclear Data File), ENDL (Evaluated Nuclear Data Library) và ACTL (the Activation

Library).

MCNP là chương trình ứng dụng phương pháp Monte Carlo để mô phỏng các quá trình vật lý

hạt nhân đối với neutron, photon, electron mang tính thống kê (các quá trình phân rã hạt nhân,

tương tác giữa các tia bức xạ với vật chất, thông lượng neutron…). MCNP sử dụng các thư viện dữ

liệu của các quá trình hạt nhân, các quy luật phân bố thống kê, số ngẫu nhiên ghi lại các sự kiện của

một hạt trong suốt quá trình kể từ khi phát ra từ nguồn đến hết thời gian sống của nó. Chương trình

này là công cụ mô phỏng được thiết lập rất đáng tin cậy, cho phép người sử dụng xây dựng các dạng

hình học phức tạp và mô phỏng dựa trên các thư viện hạt nhân. Phương pháp Monte Carlo có thể

được sử dụng để sao chụp một cách lý thuyết một quá trình thống kê (chẳng hạn như tương tác

gamma qua vật chất) và rất hữu hiệu trong các bài toán phức tạp. Chương trình MCNP điều khiển

các quá trình này bằng cách gieo số theo quy luật thống kê cho trước và mô phỏng được thực hiện

trên máy tính vì số lần thử cần thiết thường rất lớn.

MCNP từ khi ra đời cho đến nay có rất nhiều phiên bản, mỗi phiên bản kế tiếp đều được cập

nhật thêm các tính năng mới :

- Phiên bản MCNP3, 3A và 3B ra đời vào thập niên 1980 tại Los Alamos.

- Năm 1993, MCNP4A ra đời.

- Năm 1997, MCNP4B xuất hiện với việc cập nhật thêm các tính năng về photon.

- Năm 2000, phiên bản 4C ra đời kèm thêm các tính năng về electron.

- Năm 2002, phiên bản MCNPX 2.4.0 ra đời.

- Năm 2003, MCNP5 xuất hiện với các mức năng lượng, chủng loại hạt được mở rộng. Đây

là phiên bản mới nhất của MCNP.

1.2.2.2. Chương trình MCNP5

• Giới thiệu

Chương trình MCNP5, cũng giống như các phiên bản khác, sử dụng việc gieo số ngẫu nhiên

tuân theo các quy luật phân bố, ghi lại quá trình sống của một hạt khi nó được phát ra từ nguồn.

Chương trình có nhiều ứng dụng như thiết kế lò phản ứng, an toàn tới hạn, che chắn và bảo vệ, vật

lý y học … Trong phạm vi của luận văn chỉ giới hạn sử dụng chương trình này để mô phỏng phổ

gamma và tính toán hiệu suất ghi của detector.

• Cấu trúc một input file của MCNP5

Phần input file của chương trình MCNP5 được xác định như sau:

Tiêu đề và thông tin về input file (nếu cần)

Cell Cards

Surface Cards

Data Cards

(Mode Cards, Material Cards, Source Cards, Tally Cards, . . .)

• Hình học trong MCNP5

Hình học của MCNP5 thể hiện là hình học có cấu hình ba chiều tùy ý. MCNP5 xử lý các

hình học trong hệ tọa độ Descartes. MCNP5 có một chương trình dựng sẵn để kiểm tra lỗi của dữ

liệu đầu vào, thêm vào đó khả năng vẽ hình học của MCNP5 cũng giúp người dùng kiểm tra các lỗi

hình học. Sử dụng các mặt biên được xác định trên các cell card và surface card, MCNP5 theo dõi

sự chuyển động của các hạt qua các hình học, tính toán các chỗ giao nhau của các quỹ đạo vết với

các mặt biên và tìm khoảng cách dương nhỏ nhất của các chỗ giao. Nếu khoảng cách tới lần va

chạm kế tiếp lớn hơn khoảng cách nhỏ nhất, hạt sẽ rời khỏi cell đang ở. Sau đó, tại điểm giao thu

được trên bề mặt, MCNP5 sẽ xác định cell tiếp theo mà hạt sẽ vào bằng cách kiểm tra giá trị của

điểm giao (âm hoặc dương) đối với mỗi mặt được liệt kê trong cell. Dựa vào kết quả đó, MCNP5

tìm được cell đúng ở phía bên kia và tiếp tục quá trình vận chuyển. Hình học trong MCNP5 được

thể hiện qua các cell card và surface card.

a. Cell Card:

Cell là một vùng không gian được hình thành bởi các mặt biên (được định nghĩa trong phần

surface card). Nó được hình thành bằng cách thực hiện các toán tử giao, hội và bù các vùng không

gian tạo bởi các mặt. Mỗi mặt chia không gian thành hai vùng với các giá trị dương và âm tương

ứng. Khi một cell được xác định, một vấn đề quan trọng là xác định được giá trị của tất cả những

điểm nằm trong cell tương ứng với một mặt biên. Giả sử rằng s = f(x,y,z) = 0 là phương trình của

một mặt trong bài toán. Đối với một điểm (x,y,z) mà có s = 0 thì điểm đó ở trên mặt, nếu s âm điểm

đó được gọi là ở bên trong mặt và được gán dấu âm. Ngược lại, nếu s dương, điểm ở bên ngoài mặt

thì được gán dấu dương. Cell được xác định bởi cell card. Mỗi cell được diễn tả bởi số cell (cell

number), số vật chất (material number), mật độ vật chất (material density), một dãy các mặt có dấu

(âm hoặc dương) kết hợp nhau thông qua các toán tử giao, hội, bù để tạo thành cell.

b. Surface Card:

Surface card được xác định bằng cách cung cấp các hệ số của các phương trình mặt giải tích

hay các thông tin về các điểm đã biết trên mặt. MCNP5 cũng cung cấp các các dạng mặt cơ bản

chẳng hạn như mặt phẳng, mặt cầu, mặt trụ, ... (có tất cả gần 30 loại mặt cơ bản) có thể được kết

hợp với nhau thông qua các toán tử giao, hội và bù. Có hai cách để xác định các thông số mặt trong

MCNP5:

– Cung cấp các hệ số cần thiết thỏa mãn phương trình mặt.

Ví dụ: P A B C D có nghĩa là xây dựng mặt phẳng Ax + By + Cz = D.

– Xác định các điểm hình học đã biết trên các mặt đối xứng quay trên một trục tọa độ.

• Dữ liệu hạt nhân

Các bảng dữ liệu hạt nhân là những phần không thể thiếu được trong chương trình MCNP5.

Ngoài việc sử dụng các bảng dữ liệu có sẵn trong MCNP5, người dùng còn có thể sử dụng các dữ

liệu được tái tạo từ các dữ liệu gốc bên ngoài thông qua một chương trình chuyển đổi chẳng hạn như

NJOY hay là các dữ liệu mới được đưa vào trong MCNP5 bởi chính bản thân người dùng. Có tất cả

9 loại dữ liệu hạt nhân trong MCNP5:

– Tương tác neutron có năng lượng liên tục

– Tương tác neutron phản ứng rời rạc

– Tương tác quang nguyên tử năng lượng liên tục

– Tương tác quang hạt nhân năng lượng liên tục

– Các tiết diện để tính liều cho neutron

– Neutron S(a,b) nhiệt

– Tương tác neutron, cặp neutron/photon, các hạt tích điện giả neutron

– Tương tác photon

– Tương tác electron

Các dữ liệu hạt nhân được đưa vào trong MCNP5 qua phần khai báo material card.

• Mô tả nguồn

MCNP5 cho phép người dùng mô tả nguồn ở các dạng khác nhau thông qua các thông số

nguồn như năng lượng, thời gian, vị trí và hướng phát nguồn hay các thông số hình học khác như

cell hoặc surface. Bên cạnh việc mô tả nguồn theo phân bố xác suất, người dùng còn có thể sử dụng

các hàm dựng sẵn để mô tả nguồn. Các hàm này bao gồm các hàm giải tích cho các phổ năng lượng

phân hạch và nhiệt hạch chẳng hạn như các phổ Watt, Maxwell và các phổ dạng Gauss (dạng theo

thời gian, dạng đẳng hướng, cosin và dọc theo một hướng nhất định).

• Tally

Trong MCNP5 có nhiều loại tally tính toán khác nhau. Người sử dụng có thể dùng các tally

(đánh giá) khác nhau liên quan đến dòng hạt, thông lượng hạt, năng lượng để lại. Các tally trong

MCNP5 đã được chuẩn hoá trên một hạt phát ra, ngoại trừ một vài trường hợp đối với nguồn tới

hạn. Các loại tally thông dụng của MCNP5 được cho trong bảng 1.

Bảng 1: Các kiểu đánh giá (tally) trong MCNP5

Kí hiệu Tally Loại hạt

F1 Cường độ dòng qua bề mặt N, P, E

F2 Thông lượng trung bình qua một bề mặt N, P, E

F4 Thông lượng trung bình qua một cell N, P, E

F5 Thông lượng tại một điểm hay vòng N, P

F6 Năng lượng để lại trung bình qua một cell N, P

F7 Năng lượng phân hạch trung bình để lại trong một cell N

F8 Sự phân bố độ cao xung trong một cell P, E

• Output Ngoài các thông tin về kết quả, output file của MCNP5 còn có các bảng chứa các thông tin

tóm tắt cần thiết cho người sử dụng để biết rõ thêm về quá trình chạy mô phỏng của MCNP5. Các

thông tin này làm sáng tỏ vấn đề vật lí của bài toán và sự thích ứng của mô phỏng Monte Carlo. Nếu

có xảy ra sai sót trong khi chạy chương trình, MCNP5 sẽ in chi tiết cảnh báo trong phần output để

người sử dụng có thể tìm và loại bỏ. Các kết quả tính toán được in ra cùng với độ lệch chuẩn. Ngoài

ra, đi kèm với các kết quả còn là một bảng phân tích chi tiết để xác định độ tin cậy của các kết quả

này.

1.3. NGÔN NGỮ LẬP TRÌNH FORTRAN

1.3.1. Giới thiệu

Fortran là ngôn ngữ lập trình bậc cao lâu đời và được dùng phổ biến hiện nay trên thế giới,

nhất là trong lĩnh vực khoa học tự nhiên và khoa học kỹ thuật. Đây là ngôn ngữ lập trình bậc cao

xuất hiện ngay buổi đầu của kỹ thuật máy tính và phát triển cho đến ngày nay, do đó Fortran có một

kho dữ liệu là các chương trình chuẩn – chương trình mẫu phong phú nhất trong số các ngôn ngữ

lập trình bậc cao. Để tăng cường các tính năng mạnh hơn, phù hợp cho việc chuyển đổi từ hệ máy

tính này sang hệ máy tính khác, Fortran đã nhiều lần được chuẩn hóa và cải tiến. Phiên bản Fortran

đầu tiên được dùng vào năm 1957 và là ngôn ngữ lập trình bậc cao đầu tiên. Lúc đầu ngôn ngữ này

được dùng cho các tính toán khoa học nhưng theo năm tháng Fortran trở thành ngôn ngữ dùng cho

nhiều mục đích lập trình khác nhau. Tiếp theo phiên bản đầu tiên ra đời vào năm 1957, các phiên

bản khác đã lần lượt xuất hiện như Fortran II (1958), Fortran III (1958), Fortran IV (1961), Fortran

66 (1966), Fortran 77 (1977), Fortran 90 (1990), Fortran 95 (1995) và phiên bản cuối cùng gần đây

nhất là Fortran 2003 được hoàn thiện vào năm 2003 và chính thức xuất bản để ứng dụng vào năm

2004. Tuy nhiên các phiên bản được nhiều người sử dụng nhất hiện nay lại là Fortran 77 và Fortran

90. Các phiên bản mới của Fortran được phát triển trên cơ sở các phiên bản cũ chủ yếu bằng cách

thêm các tính năng mới và tiện ích mới mạnh hơn nhưng các câu lệnh cơ bản vẫn không đổi. Do đó,

việc sử dụng chương trình viết trên các phiên bản khác nhau là không gặp trở ngại.

1.3.2. Giải bài toán bằng ngôn ngữ lập trình Fortran

Để giải một bài toán bằng lập trình với ngôn ngữ Fortran, cần thực hiện theo trình tự các

bước sau:

1) Phân tích bài toán, xác định thuật giải, các bước thực hiện và trình tự thực hiện các bước.

Đây là bước hết sức quan trọng, vì nó quyết định sự đúng đắn về mặt lôgic của việc giải bài toán.

Do đó, nói chung cần phải lập một dàn bài cụ thể và biểu diễn nó qua các sơ đồ (thường gọi là sơ đồ

khối).

2) Soạn thảo mã nguồn của chương trình (chương trình nguồn, hay lời chương trình), tức là

ngôn ngữ hoá các thuật giải, theo đúng trình tự đã lập và lưu vào một (hoặc một số) file với phần

mở rộng là *.f90 (hoặc *.f, *.for, ngầm định đối với Fortran 77).

3) Tiến hành biên dịch chương trình. Ở bước này nếu chương trình vẫn còn lỗi cú pháp thì

quay lại bước 2 để chỉnh sửa rồi tiếp tục biên dịch lại chương trình. Quá trình cứ tiếp diễn cho đến

khi trình biên dịch tạo được file đích (Ojective file) và thực hiện liên kết (link) để nhận được file

thực hiện (executable file).

4) Chạy chương trình (tức chạy file thực hiện) để nhận được kết quả. Sau khi nhận được kết

quả tính toán, cần phân tích, xem xét tính hợp lý, đúng đắn của nó. Nếu kết quả không phù hợp cần

phải xem xét lại bước 1 và bước 2.

1.3.3. Cấu trúc chung của một chương trình Fortran

Cấu trúc chung của một chương trình Fortran đơn giản như sau:

PROGRAM TenChuongTrinh

Cac cau lenh khai bao

Cac cau lenh thuc hien

END PROGRAM TenChuongTrinh

TenChuongTrinh là tên của chương trình, thường được đặt một cách tùy ý sao cho mang tính

gợi nhớ, rằng chương trình sẽ giải quyết vấn đề gì.

Cac cau lenh khai bao là những câu lệnh khai báo biến, hằng,... và kiểu dữ liệu tương ứng

của chúng để trình biên dịch cấp phát bộ nhớ, phân luồng xử lý. Các câu lệnh khai báo nằm trong

nhóm các câu lệnh không thực hiện (nonexecutable statement).

Cac cau lenh thuc hien là những câu lệnh xác định qui tắc và trình tự thực hiện tính toán, xử

lý để đạt được kết quả.

PROGRAM là câu lệnh tùy chọn, không nhất thiết phải có. Nếu không được viết ra, chúng

có thể được chương trình dịch tự thêm vào.

Trong cấu trúc trên, các mục (nếu có) bắt buộc phải xuất hiện theo trình tự như đã mô tả. Có

nghĩa là sau câu lệnh mô tả tên chương trình sẽ là các câu lệnh khai báo, tiếp theo là các câu lệnh

thực hiện. Câu lệnh END phải đặt ở cuối chương trình.

Chỉ có một câu lệnh bắt buộc trong chương trình Fortran là END. Câu lệnh này báo cho

chương trình dịch rằng không còn câu lệnh nào hơn nữa để dịch. Có thể bỏ qua TenChuongTrinh

trong câu lệnh END, nhưng nếu có TenChuongTrinh thì từ khoá PROGRAM là bắt buộc.

Những dòng bắt đầu với dấu chấm than (Fortran 90) hoặc chữ “c” (Fortran 77) là lời giải

thích, làm rõ cho người đọc chương trình, và không ảnh hưởng gì tới chương trình dịch. Các biến

trong chương trình có thể có các kiểu (type) khác nhau. Các dòng trống (nếu có) trong chương trình

được xem như những câu lệnh không thực hiện (non-executable), tức là không có tác động nào được

thực hiện, có thể chèn thêm vào để cho chương trình được sáng sủa, không rối mắt.

1.3.4. Cấu trúc chung của một câu lệnh Fortran

Dạng câu lệnh cơ bản của mọi chương trình Fortran 90 có thể gồm từ 0 đến 132 ký tự (câu

lệnh có thể là trống rỗng; câu lệnh trống rỗng làm cho chương trình dễ đọc hơn bởi sự phân cách

lôgic giữa các đoạn). Đối với phiên bản Fortran 77 và các phiên bản trước đó, nội dung các câu lệnh

phải bắt đầu từ cột thứ 7 và kéo dài tối đa đến cột thứ 72. Nếu câu lệnh có nội dung dài hơn, nó sẽ

được ngắt xuống dòng dưới, và ở dòng nối tiếp này phải có một ký tự bất kỳ (khác dấu cách) xuất

hiện ở cột thứ 6. Fortran 90 không có sự hạn chế đó.

Một câu lệnh cũng có thể có nhãn. Nhãn là một số nguyên dương trong khoảng 1 đến 99999.

Nhãn (nếu có) phải là duy nhất trong một chương trình và phải đặt ở đầu câu lệnh, phân cách với

nội dung câu lệnh bởi ít nhất một dấu cách. Đối với Fortran 77 và các phiên bản trước, nhãn được

ghi vào các cột 1 đến 5.

Tất cả các câu lệnh, trừ câu lệnh gán, đều bắt đầu bằng các từ khoá (keyword). Ví dụ một số

từ khoá thường gặp như PROGRAM, OPEN, PRINT, READ, END...

CHƯƠNG 2: MÔ HÌNH HÓA HỆ PHỔ KẾ GAMMA BẰNG CHƯƠNG

TRÌNH MCNP5

Ngày nay khó có thể tìm thấy một lĩnh vực nào trong khoa học hạt nhân mà ở đấy đo phổ

gamma với việc sử dụng detector bán dẫn lại không có mặt. Đo phổ gamma trở thành công nghệ

phát triển. Trong nhiều lĩnh vực của khoa học hạt nhân ứng dụng, những tính chất của bức xạ được

sử dụng để đo nồng độ phóng xạ, chẳng hạn như xác định hàm lượng của các hạt nhân phóng xạ

phát gamma trong các mẫu môi trường. Những detector bán dẫn đã được đặt ở những vị trí chính

yếu trong các phòng thí nghiệm phân tích phóng xạ. Việc sử dụng các detector bán dẫn đã giúp tạo

nên các kết quả chính xác hơn cho việc ghi nhận các bức xạ gamma của detector với các năng lượng

khác nhau.

Có nhiều loại detector khác nhau nhưng tất cả đều dựa trên cùng một nguyên tắc là chuyển

một phần hay toàn bộ năng lượng bức xạ trong detector thành xung điện.

2.1. MÔ HÌNH HÓA HỆ PHỔ KẾ GAMMA

Để xây dựng một input (bộ số liệu đầu vào của chương trình mô phỏng MCNP5), trước hết

cần mô hình hóa hệ phổ kế gamma. Các hệ phổ kế được sử dụng trong các phòng thí nghiệm của

các trường đại học và các trung tâm nghiên cứu hiện nay về cơ bản có cấu trúc giống nhau, chỉ khác

nhau về thông số từng bộ phận. Hệ phổ kế được mô tả trong luận văn này là hệ phổ kế của Trung

tâm Hạt nhân TPHCM.

Để mô hình hoá hệ phổ kế gamma bằng chương trình MCNP5, cần phải xây dựng một cách

chính xác bộ số liệu đầu vào. Bộ số liệu này chứa đựng những thông tin về cấu trúc hình học, thành

phần vật liệu của detector, buồng chì và các nguồn phóng xạ. Đầu tiên việc mô hình hoá hệ phổ kế

được thực hiện dựa trên các chi tiết cấu hình của hệ đo, vật liệu, các thông số về mật độ, thành phần

hóa học, nồng độ từng nguyên tố tham gia trong chất cấu thành vật liệu nền tương ứng, các đặc

trưng của nguồn phóng xạ, loại phân bố năng lượng, xác suất phát hạt, loại hạt gây tương tác trên

detector, kiến thức về quá trình tương tác của các hạt quan tâm, các ảnh hưởng liên quan, loại đánh

giá cần xác định,... Việc hiểu biết tỉ mỉ mô hình giúp người sử dụng xây dựng input file chính xác.

Trong đó thông tin về cấu trúc hình học và thành phần vật liệu của buồng chì có được bằng cách

khảo sát đo đạc trực tiếp, thông tin của detector và các nguồn phóng xạ có được bằng cách dựa vào

số liệu do nhà sản xuất cung cấp. Bộ số liệu đầu vào này sẽ được đưa vào trong một input chuẩn của

chương trình MCNP5. Do đó input cần phải được chuẩn bị cẩn thận và phải đáp ứng chính xác các

yêu cầu, các khuôn mẫu của chương trình MCNP5, khi đó chương trình MCNP5 sẽ tái tạo lại mô

hình chính xác nhất trên máy tính về hệ phổ kế gamma trong phòng thí nghiệm. Sau đó dựa vào

những dữ liệu về tính chất hạt nhân và các quy luật tương tác hạt nhân từ thư viện dữ liệu nguồn của

chương trình, MCNP5 sẽ cho kết quả phổ gamma dựa trên mô phỏng Monte Carlo.

Khối xử lý Khối xử Tiền Bộ Detector HPGe biên độ đa lý và lưu kh ế h kh ế h

Nguồn

cung cấp

ế

Hình 2.1: Sơ đồ khối hệ phổ kế gamma dùng detector HPGe

Hệ phổ kế gamma tại Trung tâm Hạt nhân TPHCM bao gồm các phần chính như sau: buồng

chì, detector HPGe GC1518, nguồn cung cấp cao thế, tiền khuếch đại nhạy điện tích, khuếch đại,

khối phân tích biên độ đa kênh, khối xử lý và lưu trữ số liệu [11]. Ảnh chụp hệ phổ kế gamma

phông thấp hiện đang hoạt động được trình bày trong Phụ lục 1.

2.1.1. Cấu trúc của buồng chì

Trong quá trình đo đạc, các đồng vị phóng xạ tự nhiên và nhân tạo phân bố xung quanh

detector sẽ tạo ra phông phóng xạ làm ảnh hưởng đến kết quả phân tích phổ gamma đo được, để

P (tương đương 0,88 cm) cường độ chùm photon có năng lượng 1000

2 với tấm chì có bề dày 10 g/cmP

2 keV đã giảm bớt một nửa, tăng bề dày của tấm chì lên 10 lần, bề dày 100 g/cmP

P (tương đương 8,8

giảm bớt phổ phông này cần phải thiết kế vật liệu che chắn thích hợp. Chì là vật liệu được chọn vì

cm) thì cường độ chùm photon có năng lượng 1000 keV giảm đi 1000 lần. Do đó tấm chì có bề dày

210

10 cm thường được sử dụng để che chắn các bức xạ phông. Tuy nhiên, cần lưu ý là chì được sử

PPb (Chu kỳ bán rã

238

PU.

dụng để che chắn phải là chì “cổ” nhằm giảm đáng kể hoạt độ phóng xạ của P

T½ = 21 năm) được tạo ra từ quá trình phân rã của P

Buồng chì được chế tạo tại Trung tâm Hạt nhân TPHCM đáp ứng tốt đối với các yêu cầu trên

[4]. Cấu trúc của buồng chì được trình bày trên hình.

Hình 2.2: Mặt cắt dọc của buồng chì, kích thước tính bằng mm

Hình 2.3: Mặt cắt dọc của buồng chì được mô hình hoá bằng chương trình MCNP5. Các chữ số chỉ

thị số thứ tự thẻ mặt được mô tả trong input của chương trình MCNP5.

Trên hình này detector HPGe GC1518 là một ống hình trụ bán kính 3,81 cm với chiều cao

nằm bên trong buồng chì là 8,40 cm.

Buồng chì có dạng hình trụ với bán kính ngoài 25 cm, cao 50 cm, bán kính trong 15 cm, cao

30 cm. Bề dày tấm chì ở các mặt trên, mặt dưới và mặt bên hình trụ bằng 10 cm.

Ở mặt dưới của nắp buồng chì có một lớp thiếc dày 0,3 cm và một lớp đồng dày 0,1 cm. Mặt

trên của đáy buồng chì có lót một lớp đồng dày 0,8 cm. Mặt trong của thành buồng chì có một lớp

thiếc dày 0,8 cm, một lớp paraffin dày 6,25 cm nửa dưới và 4,75 cm nửa trên, và một lớp đồng dày

0,6 cm kể từ bên ngoài vào.

210

PPb phát ra tia X có năng lượng 46,5 keV ...) nên các

Buồng chì có thể chứa các đồng vị phóng xạ phát ra các tia X đặc trưng hoặc tia bêta bị hãm

và phát ra bức xạ bremsstrahlung (chẳng hạn P

lớp thiếc, đồng được lót thêm để hấp thụ các tia phóng xạ này.

Ngoài ra trong buồng chì còn lót thêm một lớp paraffin dùng để hấp thụ các neutron có

nguồn gốc từ vũ trụ hoặc do sự phá vỡ các hạt nhân nặng (chẳng hạn 1 kg chì có thể tạo ra 0,11

neutron/phút) là nguyên nhân gây ra phản ứng (neutron, gamma).

Lớp thiếc có bề dày 1 mm có thể hấp thụ đến 95% các tia X của chì và lớp đồng có bề dày

1,5 mm thêm vào có thể hấp thụ đến 98% các tia X của chì trong dải năng lượng từ 75 - 85 keV

[10]. Cấu trúc gồm các lớp thiếc, đồng và paraffin được lót thêm vào đã làm giảm đáng kể các tia X

trong dải năng lượng 70 - 90 keV [1].

2.1.2. Detector HPGe GC1518

Hình 2.4: Mặt cắt dọc của Detector HPGe GC1518, kích thước tính bằng mm

Hình 2.5: Mặt cắt dọc của Detector HPGe GC1518 được mô hình hoá bằng chương trình MCNP5.

Các chữ số chỉ thị số thứ tự thẻ mặt được mô tả trong input của chương trình MCNP5

Sơ đồ cấu trúc của detector HPGe GC1518 được trình bày trên hình 2.4 và hình 2.5. Đây là

detector germanium siêu tinh khiết có dạng hình trụ đồng trục. Detector có các thông số danh định

như sau:

60

PCo.

- Hiệu suất tương đối 15% so với detector nhấp nháy NaI (Tl) kích thước 7,62 cm × 7,62 cm.

60

- Độ phân giải năng lượng 1,8 keV tại vạch năng lượng 1332 keV của đồng vị P

PCo.

- Tỷ số đỉnh/Compton bằng 45:1 cũng tại vạch năng lượng 1332 keV của đồng vị P

Các thông số về cấu trúc hình học và thành phần vật liệu của detector do nhà sản xuất cung

10

cấp [11]. Phần chính của detector HPGe GC1518 là tinh thể germanium siêu tinh khiết (mức độ tạp

3 P nguyên tử/cmP

P) có đường kính ngoài 54 mm, chiều cao 32 mm, ở giữa

chất thuần vào khoảng 10P

có một hốc hình trụ đường kính 7 mm và chiều cao 17 mm. Tín hiệu lấy ra từ một điện cực bằng

P được

+ cùng với bề dày tương đương 0,35 mm. Ge được gọi là lớp germanium bất hoạt. Đó là lớp nP

đồng đặt ở trong hốc của tinh thể. Mặt trên và mặt bên tinh thể có một lớp lithium khuếch tán ngoài

+ P mm Ge. Đây là lớp pP

P được nối với cực âm của nguồn điện.

-3 dày tương đương 3.10P

nối với cực dương của nguồn điện. Mặt trong hốc của tinh thể có một lớp boron được cấy ion với bề

-3 0,1 mm và lớp dưới làm bằng mylar được kim loại hóa với bề dày 8,5.10P

P mm.

Mặt trên cùng của tinh thể có phủ hai lớp vật liệu gồm lớp trên làm bằng kapton với bề dày

Tinh thể germanium đặt trong một hộp kín bằng nhôm và ghép cách điện với que tản nhiệt

PC (77 K ) nhằm giảm tối thiểu ảnh hưởng nhiễu do dao động

0 đến bình chứa nitrogen lỏng -196 P

bằng đồng nhưng vẫn đảm bảo sự tản nhiệt tốt. Que tản nhiệt sẽ dẫn nhiệt từ tinh thể germanium

nhiệt trong tinh thể germanium và các linh kiện điện tử của tiền khuếch đại.

Hộp kín bằng nhôm có bề dày 2,7 mm ở chỗ dày nhất và 0,76 mm ở chỗ mỏng nhất để đảm

bảo tránh được sự hấp thụ các photon có năng lượng thấp và che chắn bức xạ hồng ngoại từ bên

ngoài vào tinh thể germanium.

Các điện cực cách điện với nhau bằng teflon và có một khoảng chân không ở dưới tinh thể.

Toàn bộ hộp kín này được đặt trong một vỏ bằng nhôm có đường kính 76,2 mm và bề dày 1,5 mm.

Khoảng chân không giữa mặt trên tinh thể germanium với mặt dưới của vỏ nhôm là 5 mm để tránh

các va chạm vào bề mặt tinh thể germanium khi lắp ráp detector. Detector HPGe GC1518 được đặt

trong buồng chì để giảm phông gamma từ môi trường.

2.1.3. Cấu trúc của các nguồn phóng xạ

Trong thực nghiệm nhiều loại nguồn phóng xạ do các hãng khác nhau sản xuất được sử dụng,

do đó cấu trúc hình học và thành phần vật liệu của chúng cũng rất đa dạng. Có loại là nguồn phóng

xạ chuẩn, hình học nguồn điểm thường dùng để hiệu chuẩn (calibration) hệ phổ kế gamma; có loại

là nguồn hỗn hợp, hình học nguồn rộng, trong đó chất phóng xạ được trộn đều với vật liệu nền như

bột đá vôi, bột zirconium silicate ...

Khi mô hình hóa hệ phổ kế gamma, việc mô phỏng nguồn phóng xạ cũng là một vấn đề rất

quan trọng, nguồn phóng xạ có nhiều loại khác nhau và cần được mô tả tỉ mỉ về cấu trúc hình học

cũng như thành phần vật liệu, từ đó mới có thể đảm bảo kết quả mô phỏng chính xác. Ở đây, chúng

137

tôi giới thiệu một nguồn phóng xạ chuẩn rất hay được sử dụng trong nghiên cứu phổ gamma là

PCs, đại diện cho phổ gamma đơn năng.

137

nguồn P

PCs của hãng Nuclear Services & Supplies - Rost GmbH đặt tại Trung tâm

Nguồn phóng xạ P

Hạt nhân TPHCM có dạng điểm đường kính 2 mm đặt trong một hốc hình giếng (giống như

collimator) của một giá đỡ bằng thép không gỉ dạng hình trụ được trình bày trong hình 2.7. Trên giá

đỡ này đầu đặt nguồn phóng xạ có dạng hình trụ chuẩn, đầu còn lại được vát phẳng ở hai mặt đối

137

diện nhằm thuận tiện trong việc lắp đặt nguồn để tiến hành thực nghiệm đo đạc. Sơ đồ phân rã của

PCs được trình bày trong hình 2.6.

đồng vị phóng xạ P

137

PCs

Hình 2.6: Sơ đồ phân rã của đồng vị phóng xạ P

Đường gạch ngang đậm: trạng thái cơ bản của hạt nhân mẹ hoặc hạt nhân con. Đường gạch ngang

nhạt: trạng thái kích thích của hạt nhân con. Cột số bên trái: spin và độ chẵn lẻ. Cột số bên phải:

137

năng lượng tính bằng MeV.

PCs của hãng Nuclear Servisces & Supplies

Hình 2.7: Cấu trúc của nguồn P

137

PCs

– Rost GmbH

Bảng 2: Các thông tin về nguồn P

137

P

Năng Cường Hoạt Nơi Ngày Tên Chu kỳ lượng độ phát độ sản xuất sản xuất nguồn bán rã (keV) xạ (%) (kBq)

PCs

30,04 năm 661,6 85,21 370 NSS, USA 30-11-1990

2.1.4. Mô phỏng phổ gamma của các nguồn phóng xạ

Các photon phát ra từ nguồn sẽ phân bố đều theo mọi hướng trong không gian. Chỉ có một

phần các photon đạt đến bề mặt detector, phần còn lại sẽ bị hấp thụ trong môi trường xung quanh hệ

phổ kế. Một photon kể từ khi phát ra từ nguồn cho đến lúc kết thúc có thể xảy ra các quá trình tương

tác với vật chất trên suốt quãng đường truyền qua của nó. Tuy nhiên chỉ có các quá trình tương tác

của photon với vật liệu bên trong thể tích germanium hoạt động mới đóng góp số đếm vào phổ

gamma. Dựa trên cơ sở các đặc điểm và chuẩn mực của chương trình MCNP5, bộ số liệu đầu vào về

cấu trúc hình học và thành phần vật liệu của buồng chì, detector và nguồn phóng xạ được đưa vào

input của chương trình MCNP5 để mô hình hoá hệ phổ kế gamma và mô phỏng phổ gamma của các

nguồn phóng xạ sao cho thời gian tính toán càng ngắn càng tốt nhưng vẫn phải đảm bảo độ tin cậy

của phổ gamma mô phỏng.

2.2. INPUT CỦA CHƯƠNG TRÌNH MCNP5

Mỗi tính toán MCNP5 cụ thể cần được cung cấp tập các số liệu đầu vào chứa đựng thông tin

liên quan đến thư viện các tiết diện và mô tả hình học vật lý của nguồn, detector và các vật liệu khác

cũng như năng lượng của gamma, các khe năng lượng (tương ứng với kênh trong phổ được đo) mà

trong đó các sự kiện được đánh giá đối với năng lượng bị mất trong thể tích detector và số các

photon được phát ra.

MCNP5 theo dõi mỗi photon khi nó truyền qua không gian và tương tác với các nguyên tử

trong các vật liệu khác nhau có mặt ở đó. Các electron và các photon thứ cấp được tạo ra trong

những tương tác này cũng được theo dõi cho đến khi toàn bộ năng lượng của chúng bị mất trong các

vật liệu khác nhau hoặc chúng thoát ra khỏi không gian vật lý được bao hàm trong mô hình.

Đối với các tương tác trong thể tích detector, MCNP5 tạo ra đánh giá số các sự kiện trong

mỗi khe năng lượng tức là nó cung cấp phổ mất mát năng lượng. Vì hệ đo không đo trực tiếp năng

lượng để lại trong detector nên phổ mô phỏng sẽ khác ở một chừng mực nào đó với phổ đo, thậm

chí nếu phổ mô phỏng không có bất kỳ sự xấp xỉ hay sai số. Đối với detector bán dẫn Ge, kiểu

detector có đặc trưng rất tuyến tính (tức là biên độ của tín hiệu từ detector tỷ lệ với năng lượng để

lại trong nó) và độ phân giải năng lượng rất tốt (tức các đỉnh quan sát được đều rất hẹp) nên sự khác

nhau này thường là rất nhỏ.

Như đã trình bày ở trên hệ phổ kế gamma gồm buồng chì, detector, nguồn phóng xạ và hệ

thống điện tử rất phức tạp. Tuy nhiên khi tiến hành mô hình hoá hệ phổ kế thì có thể bỏ qua những

phần không gian đóng góp không đáng kể vào phổ gamma mô phỏng [6], do đó chỉ cần mô tả cấu

trúc hình học và thành phần vật liệu của buồng chì, detector và nguồn phóng xạ.

Trong quá trình mô phỏng mode p được sử dụng bởi vì nguyên tố germanium có bậc số

nguyên tử Z lớn cho nên sự khác nhau giữa mode p và mode p e là không đáng kể [2]. Mặt khác mô

hình chi tiết về tương tác của photon với vật chất cũng được áp dụng, trong mô hình này ngoài việc

tính toán đối với các quá trình tương tác quan trọng như hấp thụ quang điện, tán xạ Compton (tán xạ

không kết hợp), hiệu ứng tạo cặp còn phải tính toán đối với quá trình tán xạ Thomson (tán xạ kết

hợp) và quá trình phát huỳnh quang xảy ra theo sau quá trình hấp thụ quang điện. Đối với mode p

quá trình tương tác của electron với vật chất được mô phỏng theo mô hình gần đúng TTB (thick

target bremsstrahlung) của chương trình MCNP5. Mô hình gần đúng TTB giả thiết rằng electron

được tạo thành di chuyển cùng hướng với photon tới và phát ra bức xạ bremsstrahlung ngay tức thì.

Khi photon đi xuyên qua vùng nghèo thì sẽ tạo ra các cặp hạt mang điện và được tập hợp về hai điện

cực. Thông qua bộ tiền khuếch đại nhạy điện tích, điện tích của các hạt mang điện chuyển đổi thành

xung điện áp. Xung điện áp tỉ lệ với phần năng lượng của photon được giữ lại trong detector. Khi đó

phổ phân bố độ cao xung hay còn gọi là phổ gamma mô phỏng được lấy ra bằng thẻ truy xuất kết

quả F8 của chương trình MCNP5. Khi được truy xuất bằng thẻ F8, kết quả phân bố độ cao xung

được tính bằng số đếm đối với năng lượng (chuẩn theo số quá trình photon phát ra từ nguồn tại năng

lượng đó). Ngoài ra do ảnh hưởng của ba hiệu ứng là sự giãn rộng thống kê số lượng các hạt mang

điện, hiệu suất tập hợp điện tích và đóng góp của các nhiễu điện tử [15] làm cho các quang đỉnh của

phổ gamma thực nghiệm có dạng Gauss. Do đó trong quá trình mô phỏng còn sử dụng lựa chọn

GEB (gaussian energy broadening) của thẻ FT8 đi kèm với thẻ kết quả phân bố độ cao xung F8. Với

lựa chọn GEB phổ gamma mô phỏng phù hợp tốt hơn với phổ gamma thực nghiệm. Trong phổ

gamma mô phỏng quang đỉnh được mở rộng bằng cách lấy mẫu ngẫu nhiên theo dạng hàm Gauss.

Input chuẩn được chia làm ba phần là định nghĩa ô, định nghĩa mặt và định nghĩa vật liệu.

Chúng được ngăn cách nhau bằng các dòng trống. Định nghĩa ô dựa trên các mặt biên được liên kết

lại với nhau tạo thành và được lấp đầy vật chất đồng nhất tương ứng. Định nghĩa mặt là các dạng

toàn phương liên kết tạo thành các ô. Trong định nghĩa dữ liệu cần phải khai báo nguồn, vật liệu cấu

tạo các ô, loại đánh giá cần tính toán, số hạt gieo, độ quan trọng của các ô. Một input điển hình của

chương trình MCNP5 được xây dựng để mô phỏng phổ gamma của các nguồn phóng xạ được trình

bày trong phụ lục 2. Trong đó:

Dòng 1 và 2 là các dòng tiêu đề và dòng thông báo bắt đầu khai báo thẻ ô tương ứng.

Các dòng từ 3 đến 28 khai báo các thẻ ô (cell card).

Dòng 29 là dòng thông báo bắt đầu khai báo thẻ mặt.

Các dòng từ 30 đến 70 khai báo các thẻ mặt (surface card).

Dòng 71 thông báo bắt đầu khai báo thẻ dữ liệu.

Các dòng còn lại khai báo các thẻ dữ liệu.

Cụ thể là dòng 72 mô tả mode p được sử dụng. Với mode p, quá trình vận chuyển của

electron được tính toán theo mô hình gần đúng TTB (thick target bremsstrahlung).

Dòng 95 và 88 mô tả thẻ truy xuất kết quả phân bố độ cao xung theo năng lượng F8 và thẻ xử

lý đặc biệt FT8 với lựa chọn GEB (gaussian energy broadening).

Dòng 99 mô tả kỹ thuật cắt năng lượng.

Dòng 101 và 102 mô tả điều kiện kết thúc quá trình mô phỏng gồm số photon phát ra từ

nguồn và thời gian tính toán tương ứng. Nếu trong input có các dòng bắt đầu bằng kí tự "c" thì các

dòng đó không có ý nghĩa hoặc tạm thời bỏ đi, khi chạy MCNP5 sẽ không xử lý các dòng này.

CHƯƠNG 3: XÂY DỰNG CHƯƠNG TRÌNH TÍNH TOÁN

3.1. TỔNG QUÁT VỀ CHƯƠNG TRÌNH AUT-MCNP5

Trong quá trình mô phỏng phổ gamma bằng chương trình MCNP5, người sử dụng chương

trình thường phải tiến hành những công đoạn rời rạc: tạo lập input, dùng MCNP5 để chạy input, sau

khi chương trình xuất ra file output thì người sử dụng sẽ tiếp tục quá trình lọc, tách những dữ liệu

cần thiết trong file output, vẽ phổ gamma, chuẩn phổ mô phỏng và phổ thực nghiệm để có thể so

sánh phổ gamma mô phỏng với phổ thu được từ kết quả thực nghiệm, tính toán các đại lượng liên

quan khác như diện tích đỉnh phổ, hiệu suất… Các công đoạn sau khi có file output thường được

thực hiện thủ công, vì vậy gây ra một số khó khăn cho người nghiên cứu, đặc biệt là trong những

trường hợp mà người nghiên cứu phải tiến hành chạy MCNP5 và xử lý một lượng lớn các file input

(ví dụ tính toán hiệu suất đối với các dạng hình học nguồn và mẫu khác nhau, khảo sát hiệu suất

theo năng lượng, theo khoảng cách, hiệu chỉnh trùng phùng tổng đối với gamma phân rã nhiều tầng,

hiệu chỉnh tự hấp thụ đối với hình học nguồn và mẫu thể tích,…). Chính vì vậy, chương trình AUT-

MCNP5 được xây dựng với mục đích giúp người sử dụng có thể sử dụng chương trình MCNP5 một

cách thuận tiện hơn trong bài toán mô phỏng phổ gamma, rút ngắn thời gian và giảm bớt khối lượng

công việc. Hơn nữa, việc tự động hóa quá trình tính toán còn giúp người sử dụng giảm thiểu sai sót,

đảm bảo tính chính xác của kết quả.

Chương trình AUT-MCNP5 được xây dựng trên cơ sở ngôn ngữ lập trình Fortran. Chương

trình tích hợp nhiều modun, mỗi modun của chương trình giải quyết một vấn đề khác nhau trong

quá trình tiến hành mô phỏng và xử lý phổ gamma mô phỏng. Mỗi modun đó ứng với một bước của

quá trình mô phỏng và xử lý thực tế, tuy nhiên, thời gian xử lý, tính toán được rút gọn đi rất nhiều,

USơ đồ khối chương trình

số liệu tính toán cũng đảm bảo độ chính xác cao hơn hẳn so với khi xử lý thủ công.

Bắt đầu

Nhập lựa chọn modun

Sai Lựa

h

Đúng

Modun tính toán

Hiển thị kết quả

Sai Chọn kết thúc

Chươ t ì h

Đúng

Kết thúc

Hình 3.1: Sơ đồ khối chương trình AUT-MCNP5

Dựa trên quy trình chạy MCNP5 và xử lý output thường gặp, chúng tôi xây dựng chương

trình gồm bảy lựa chọn modun tính toán và một sự lựa chọn để thoát khỏi chương trình. Mỗi modun

tính toán có một nhiệm vụ riêng.

1. Chạy MCNP5 tự động

2. Trích xuất dữ liệu từ file output

3.Trích xuất dữ liệu từ file thực nghiệm

4. Chuẩn phổ Modun tính toán

5. Tính diện tích đỉnh phổ và hiệu suất

6. Vẽ hình hệ phổ kế

7. Hình ảnh phổ gamma

Hình 3.2: Các modun tính toán của chương trình AUT-MCNP5

Hình 3.3: Cửa sổ chương trình AUT-MCNP5

Muốn sử dụng modun nào thì chọn số tương ứng, ví dụ muốn chạy MCNP5 thì nhập 1 rồi

nhấn Enter.

3.2. CÁC MODUN TÍNH TOÁN CỦA CHƯƠNG TRÌNH AUT-MCNP5

3.2.1. Modun 1 - Chạy MCNP5 tự động

Sau khi đã thiết lập các file input, nhiệm vụ tiếp theo là đưa các file input này vào chương

trình MCNP5 để chạy và xuất ra các file kết quả (file output). Quá trình chạy chương trình thường

được thực hiện thủ công bằng cách đưa lần lượt từng file input vào chương trình MCNP5. Việc này

sẽ gặp khó khăn nếu phải chạy nhiều hay rất nhiều file input (phải chờ chạy xong file này rồi tiếp

tục lặp lại thao tác để chạy file khác trong khi thời gian chạy một file input vào khoảng vài giờ đến

vài chục giờ). Vì lý do đó, chúng tôi viết chương trình chạy tự động MCNP5, chương trình này sẽ tự

động chạy các file input được chỉ định và xuất ra file output. Nhờ vậy, người sử dụng chỉ cần thực

hiện một số thao tác ban đầu mà không phải chờ đợi để lặp lại thao tác cho từng file input riêng lẻ.

Điều này rất tiện ích trong quá trình mô phỏng phổ gamma với hàng trăm, hàng ngàn file input được

đưa vào chương trình, việc chạy chương trình sẽ không bị ngắt quãng và người sử dụng có thể hoàn

toàn chủ động về mặt thời gian.

Vì thao tác “chạy chương trình MCNP5” được lặp lại nhiều lần nên trong modun này có tích

hợp chương trình con. Chương trình con “Chạy MCNP5” có nhiệm vụ quan trọng nhất là kết nối

với MCNP5 thông qua cửa sổ chương trình Fortran chính, từ đó các file input cần xử lý được tự

động đưa vào MCNP5 một cách tuần tự. Chương trình sẽ chạy liên tục cho đến khi xử lý hết số file

input mà người sử dụng yêu cầu.

USơ đồ khối chương trình

Bắt đầu

Nhập lựa chọn chạy bao nhiêu file

Sai 1 file

Đúng Nhập tên file chứa tên các file input

Nhập tên file input

Đọc lần lượt các file input

Kết nối với chương trình con chạy MCNP5

- Xuất ra file output

Kết thúc

Hình 3.4: Sơ đồ khối của chương trình “Chạy MCNP5 tự động”

USơ đồ khối của chương trình con

Bắt đầu

Kết nối chương trình Fortran và MCNP5

-> Chạy MCNP5

Xuất ra các file output

Thực hiện vòng lặp với các file input tiếp theo

Kết thúc

Hình 3.5: Sơ đồ khối của chương trình con “Chạy MCNP5 tự động”

Sơ đồ khối trên thể hiện thuật giải của chương trình, dựa vào thuật giải đó, chúng tôi viết

UĐoạn chương trình

thành đoạn chương trình như sau:

Program Chay tu dong MCNP5

use dflib

use portlib

use msflib

character*10 inputfile(100),tentaptin,ten

integer*2 N

print*,'***********************************'

print*,'* CHUONG TRINH CHAY TU DONG MCNP5 *'

print*,'***********************************'

print*,''

print*,'Ban muon chay MCNP5 bao nhieu file?'

read*,N

if(N.eq.1) then

print*,'Nhap ten file input'

read*,ten

call chayMCNP(ten)

else

print*,'Nhap ten tap tin chua ten cac file input'

read*,tentaptin

open(1,file=tentaptin)

do 1 i=1,N

read(1,'(a10)') inputfile(i)

1 enddo

do 2 j=1,N

call chayMCNP(inputfile(j))

2 enddo

goto 10

stop

endif

goto 10

UChương trình con

subroutine chayMCNP(ten)

use dflib

logical result2

character*10 ten

character*50 x

x='mcnp name='//ten//'xsdir=xsdir'

result2=systemqq(x)

return

end Đoạn chương trình trên nhằm thực hiện chạy MCNP5 tự động đối với file input được chỉ

định. Nếu muốn chạy tự động nhiều file chỉ cần thiết lập một file chứa tên tất cả các file input người

sử dụng muốn chạy. Thành công trong việc kết nối Fortran với MCNP5 sẽ là tiền đề để thiết lập một

chương trình liên tục từ chạy MCNP5, phân tích, xử lý output đến tính toán những đại lượng quan

trọng trong quá trình phân tích phổ gamma. Tất cả sẽ được kết nối với nhau và được xử lý tự động,

nhờ đó việc khảo sát những bài toán cần nhiều quá trình xử lý MCNP5 sẽ được tiến hành nhanh

chóng, chính xác và hiệu quả cao, thậm chí đối với những yêu cầu nghiên cứu trong đó phải tiến

UCửa sổ chương trình

hành phân tích phổ gamma hàng ngàn lần.

Hình 3.6: Cửa sổ chương trình “Chạy tự động MCNP5” – Nhập số file

Nếu muốn chạy MCNP5 một file, nhập 1 rồi nhấn Enter. Chương trình sẽ tiếp tục hỏi người

sử dụng tên của file input.

Hình 3.7: Cửa sổ chương trình “Chạy tự động MCNP5” – Nhập tên file input

Nhập tên file input muốn xử lý rồi nhấn Enter (ví dụ cs40g), tên file này không quá 7 ký tự,

file input phải được chép trong cùng folder chứa chương trình “Sử dụng hiệu quả MCNP5”. File

output tạo ra trong folder sẽ có tên “cs40go”.

Nếu muốn chạy nhiều file input, tạo một file chứa tên các file input, trong file này có thể

chứa nhiều hơn số file cần chạy. Ví dụ tạo một file có tên “fileinput” trong đó ghi tên 20 file input

u1, u2, u3,…, u20 (xuống dòng để ngăn cách tên các file).

Hình 3.8: Cửa sổ chương trình “Chạy tự động MCNP5” – Nhập tên tập tin chứa tên các file input

Trong cửa sổ trên nhập tên “fileinput” và nhấn Enter. Sau khi chạy chương trình, sẽ có ba file

output được tạo thành có tên là u1o, u2o và u3o.

3.2.2. Modun 2 - Trích xuất dữ liệu từ file output

Sau khi đã có được file output, người nghiên cứu thường chỉ cần lấy ra một số dữ liệu nhất

định. Trong vấn đề mô phỏng phổ gamma, những dữ liệu cần quan tâm gồm có số đếm tương ứng

với các kênh năng lượng (đối với MCNP5, số đếm này thể hiện dưới dạng phân bố độ cao xung vi

phân), tổng số hạt đầu vào để khảo sát, tổng số hạt để lại năng lượng khác 0 trong detector, tổng giá

USơ đồ khối chương trình

trị độ cao xung năng lượng.

Bắt đầu

Nhập lựa chọn trích xuất dữ liệu từ bao nhiêu file output

Sai

1 file

Đúng

Nhập file chứa các tên file output Nhập tên file output

Đọc lần lượt các file output

Trích xuất dữ

liệu từ file output Kết nối chương trình con trích

xuất dữ liệu từ file output

Xuất ra các file dữ liệu từ output

Kết thúc

Hình 3.9: Sơ đồ khối chương trình “Trích xuất dữ liệu từ file output”

Trong trường hợp chỉ xử lý một file output, chương trình sẽ trích xuất dữ liệu trực tiếp.

Trong trường hợp cần xử lý nhiều file output, chương trình sẽ kết nối với chương trình con để tiến

hành xử lý lần lượt từng file output.

USơ đồ khối của chương trình con

Bắt đầu

Mở file output

Xác định vị trí dòng bắt đầu đọc dữ liệu

Đọc các dòng dữ liệu được xác định, mỗi dòng lấy ra ba giá trị

tương ứng với kênh năng lượng, số đếm và sai số

Đọc dòng được xác định, lấy ra giá trị tổng số đếm ứng với tất cả

các kênh (dưới dạng phân bố độ cao xung năng lượng) có trừ đi

ầ

Xác định vị trí dòng, vị trí cột và lấy ra

tổng số hạt khảo sát

Xác định vị trí dòng, vị trí cột và lấy ra giá trị

tổng độ cao xung năng lượng

Xác định vị trí dòng, vị trí cột và lấy ra tổng số hạt

để lại năng lượng khác 0 trong detector

Thực hiện vòng lặp với các file output tiếp theo

Kết thúc

Hình 3.10: Sơ đồ khối chương trình con“Trích xuất dữ liệu từ file output”

Trong chương trình “Trích xuất dữ liệu từ file output”, chương trình chính sẽ đọc tên các file

output cần lấy dữ liệu và đưa vào chương trình con. Chương trình con có nhiệm vụ đọc từng file

output và xác định chính xác vị trí các số liệu cần lấy ra. Vì các file output có cùng định dạng nên để

xác định vị trí các dữ liệu này, chương trình sẽ xác định vị trí dòng, vị trí cột bắt đầu của dữ liệu và

UĐoạn chương trình

độ dài dữ liệu.

Program trichxuatoutput

use dflib

use portlib

use msflib

character*10 tenfile,out(100),tenoutput,filexuat

character*13 strnps

real*4 e2(100,10000),de(100,10000),re(100,10000),nat22,

+ nat(100),toteng(100),toteng22,e22(10000),

+ de22(10000),re22(10000)

integer nht(100),nps(100),T,nht22,nps22

character*1 aa2,bb2,cc2

print*,'************************************************'

print*,'*CHUONG TRINH TRICH XUAT DU LIEU TU FILE OUTPUT*’

print*,'************************************************'

print*,''

print*,'Ban muon trich xuat du lieu bao nhieu file?'

read*,T

select case(T)

case(1)

print*,'Nhap ten file output'

read*,tenoutput

print*,'Nhap ten file xuat ra'

read*,filexuat

open(15,file=tenoutput,form='formatted')

bb2='e'

cc2='y'

101 if(.not.eof(15)) read(15,'(130a1)') (dummy(l),l=1,130)

if((dummy(7).ne.bb2).or.(dummy(12).ne.cc2)) goto 101

do 1515 k=1,8194

1515 read(15,'(2en14.3,f7.4)') e22(k),de22(k),re22(k)

read(15,'(14a1,en14.3,7a1)')

+(dum(l),l=1,14),toteng22,(dum(l),l=1,7)

toteng22=toteng22-de22(1)-de22(2)

aa2='1'

2020 if(.not.eof(15)) read(15,'(130a1)') (dummy(l),l=1,130)

if(dummy(1).ne.aa2) goto 2020

write(strnps,*) (dummy(l),l=91,102)

read(strnps,'(i13)') nps22

do 2525 k=1,2

2525 read(15,'(130a1)') (dummy(l),l=1,130)

read(15,'(36a1,en12.3,54a1)')

+(dum(l),l=1,36),nat22,(dum(l),l=1,54)

do 3535 k=1,3

3535 read(15,'(130a1)') (dummy(l),l=1,130)

read(15,'(36a1,i12)') (dum(l),l=1,36),nht22

close(15)

open(99,file=filexuat)

write(99,7373)

7373 format('dien tich dinh nang luong toan phan:')

write(99,'(f12.7)') toteng22

write(99,7474)

7474 format('tong phan bo do cao xung nang luong:')

write(99,'(f12.7)') nat22

write(99,7575)

7575 format('tong so hat khao sat ban dau:')

write(99,'(i12)') nps22

write(99,7676)

7676 format('tong so hat bo nang luong khac 0 trong detector:')

write(99,'(i12)') nht22

write(99,787)

787 format('Kenh nang luong (MeV) Do cao xung Sai so')

do 7272 k=1,8194

7272 write(99,6868) e22(k),de22(k),re22(k)

6868 format(f12.7,5x,en14.3,7x,f7.4)

close(99)

goto 10

case(2:)

print*,'Nhap file chua ten cac file output'

read*,tenfile

open(8,file=tenfile)

do 66 j=1,T

read(8,'(a10)') out(j)

66 enddo

do 69 j=1,T

call extract_cal(j,out,e2,de,re,toteng,nat,nht,nps)

69 enddo

do 71 j=1,T

write(ordnum,'(i0)') j

open(9,file='extract'//ordnum,status='unknown')

write(9,73)

73 format('dien tich dinh nang luong toan phan:')

write(9,'(f12.7)') toteng(j)

write(9,74)

74 format('tong phan bo do cao xung nang luong:')

write(9,'(f12.7)') nat(j)

write(9,75)

75 format('tong so hat khao sat ban dau:')

write(9,'(i12)') nps(j)

write(9,76)

76 format('tong so hat bo nang luong khac 0 trong detector:')

write(9,'(i12)') nht(j)

write(9,78)

78 format('Kenh nang luong (MeV) Do cao xung Sai so')

do 72 k=1,8194

72 write(9,68) e2(j,k),de(j,k),re(j,k)

68 format(f12.7,5x,en14.3,7x,f7.4)

71 close(9)

end select

goto 10

stop UChương trình con

Subroutine extract_cal(j,out,e2,de,re,toteng,nat,nht,nps)

implicit real*8 (a-h,o-z)

save

character*10 out(100)

character*1 dummy(200),dum(200),aa,bb,cc

character*13 strnps

real*4 e2(100,10000),de(100,10000),re(100,10000),

+ toteng(100),nat(100)

integer nht(100),nps(100)

open(5,file=out(j),form='formatted')

bb='e'

cc='y'

10 if(.not.eof(5)) read(5,'(130a1)') (dummy(l),l=1,130)

if((dummy(7).ne.bb).or.(dummy(12).ne.cc)) goto 10

do 15 k=1,8194

15 read(5,'(2en14.3,f7.4)') e2(j,k),de(j,k),re(j,k)

read(5,'(14a1,en14.3,7a1)') (dum(l),l=1,14),

+ toteng(j),(dum(l),l=1,7)

toteng(j)=toteng(j)-de(j,1)-de(j,2)

aa='1'

20 if(.not.eof(5)) read(5,'(130a1)') (dummy(l),l=1,130)

if(dummy(1).ne.aa) goto 20

write(strnps,*) (dummy(l),l=91,102)

read(strnps,'(i13)') nps(j)

do 25 k=1,2

25 read(5,'(130a1)') (dummy(l),l=1,130)

read(5,'(36a1,en12.3,54a1)') (dum(l),l=1,36),

+ nat(j),(dum(l),l=1,54)

do 35 k=1,3

35 read(5,'(130a1)') (dummy(l),l=1,130)

read(5,'(36a1,i12)') (dum(l),l=1,36),nht(j)

close(5)

return

end

UCửa sổ chương trình

Hình 3.11: Cửa sổ chương trình “Trích xuất dữ liệu từ file output”-Nhập số file

Nhập số file output muốn lấy dữ liệu.

Hình 3.12: Cửa sổ chương trình “Trích xuất dữ liệu từ file output” - Nhập

tên file output

Nếu người sử dụng muốn trích xuất dữ liệu từ một file output, chương trình sẽ hỏi tên file

này. Trong ví dụ trên, tên file output được nhập vào là cs40go. Cần chú ý là tên file này chứa không

quá 10 ký tự.

Hình 3.13: Cửa sổ chương trình “Trích xuất dữ liệu từ file output” - Nhập

tên file xuất ra

Người sử dụng nhập tên file chứa các dữ liệu trích xuất được từ file output. Nếu muốn file

xuất ra có định dạng là file gì thì chọn phần mở rộng tương ứng. Ví dụ muốn file xuất ra có định

dạng file Excel thì nhập tên file xuất ra có phần mở rộng là “.xls” (cs.xls).

Hình 3.14: Cửa sổ chương trình “Trích xuất dữ liệu từ file output”-Nhập file chứa tên các file

output

Trong trường hợp người sử dụng muốn trích xuất dữ liệu nhiều file output thì liệt kê tên các

file output cần lấy dữ liệu vào một file, ví dụ “fileoutput”, xuống dòng để ngăn cách các tên file

output. Nhập “fileoutput” vào dấu nháy trên, khi đó chương trình sẽ trích xuất dữ liệu từ các file

output và xuất ra các file extract1, extract2, extract3,… Tên các file xuất đã được mặc định theo thứ

tự các file output được liệt kê. Người sử dụng có thể đổi tên các file xuất ra sau khi đã chạy xong

chương trình.

3.2.3. Modun 3 - Trích xuất dữ liệu từ file thực nghiệm

Để xác định độ tin cậy của kết quả mô phỏng bằng chương trình MCNP5, cần phải tiến hành

so sánh với kết quả thu được từ thực nghiệm. File kết quả ứng với mỗi chương trình hiển thị số liệu

trong thực nghiệm có định dạng cố định, do đó có thể thiết lập chương trình sắp xếp lại dữ liệu để

USơ đồ khối chương trình

tiện lợi cho việc tính toán và vẽ phổ gamma sau đó.

Nhập lựa chọn xử lý phổ thực nghiệm nào Bắt đầu

Đúng Sai

Chọn 1 (Phổ TT hạt nhân)

Nhập lựa chọn xử lý bao nhiêu file Nhập lựa chọn xử lý bao nhiêu file

Sai Sai

1 file 1 file

Đúng Nhập tên file Đúng Nhập tên file

chứa tên các file chứa tên các file Nhập tên file Nhập tên file thực nghiệm thực nghiệm Đọc lần lượt các Đọc lần lượt các file thực nghiệm file thực nghiệm

Kết nối với chương trình Kết nối với chương Xử lý file Xử lý file con xử lý phổ TT hạt nhân trình con xử lý phổ Bộ

Xuất file kết quả Kết thúc

USơ đồ khối chương trình con

Hình 3.15: Sơ đồ khối chương trình “Trích xuất dữ liệu từ file thực nghiệm”

Bắt đầu

Mở file kết quả thực nghiệm

- Bỏ các dòng không chứa dữ liệu trong file kết quả thực nghiệm

Xác định vị trí bắt đầu đọc dữ liệu

Xác định vị trí, khoảng cách các cột dữ liệu để đọc

dữ liệu (số đếm ứng với các mức năng lượng)

Định dạng dữ liệu xuất ra

Thực hiện vòng lặp với các file kết quả thực nghiệm tiếp theo

Kết thúc

Hình 3.16: Sơ đồ khối chương trình con “Trích xuất dữ liệu từ file thực nghiệm”

Các file kết quả thực nghiệm từ các chương trình khác nhau cũng như các hệ đo khác nhau sẽ

có định dạng khác nhau. Trong mỗi file, số đếm ứng với các kênh năng lượng được chia thành nhiều

cột dữ liệu với những khoảng cách nhất định giữa các cột dữ liệu đó. Trong chương trình này, chúng

tôi chia ra hai trường hợp là file thực nghiệm có định dạng khác nhau từ hệ phổ kế của Trung tâm

Hạt nhân TPHCM và từ Phòng thí nghiệm của Bộ môn Vật lý hạt nhân Đại học Khoa học tự nhiên

UĐoạn chương trình

TPHCM.

Program trichxuatthucnghiem

use dflib

use portlib

use msflib

character*1 dummy(200),aa

character*6 dumm,datfile2(100),ordnum2

dimension ib(10000),e(10000),ic(10000),ec(10000)

character*20 ten1,ten2,tentaptin3,ten4,ten5,tentaptin33,

+ datfile2(100),datfile(100)

integer N3,M,L3

character*4 datfile(100),ordnum

dimension ia(100,10000),id(100,10000)

integer nia(100,10000),nid(100,10000)

print*,'***************************************************'

print*,'CHUONG TRINH TRICH XUAT DU LIEU TU FILE THUC NGHIEM'

print*,'***************************************************'

print*,''

print*,'Ban muon xu ly pho thuc nghiem nao?'

print*,'1- Trung tam Hat nhan TP.HCM'

print*,'2- Bo mon VLHN Dai hoc KHTN'

read*,N3

if(N3.eq.1) then

print*,'Ban muon xu ly bao nhieu file?'

read*,M

if(M.eq.1) then

print*,’Nhap ten file thuc nghiem’

read *,ten1

open(3,file=ten1,form='formatted')

print*,’Nhap ten file xuat ra’

read *,ten2

open(33,file=ten2)

aa=':'

do i=1,9

read(3,'(76a1)')(dummy(j),j=1,76)

enddo

i=0

3 if (.not.eof(3)) then

read(3,'(6a1)') (dummy(j),j=1,6)

if(dummy(6).ne.aa) goto 3

backspace (3)

read(3,33) dumm,(ib(i+l), l=1,5)

i=i+5

max=i

goto 3

goto 10

else

do i=1,max

e(i) = .005616 + 0.000236*(i-1)

write(33,333) e(i),ib(i)

enddo

goto 10

endif

3030 format(a1)

33 format(a6,5i13)

333 format(f12.6,5i13)

else

print*,'Nhap file chua ten cac file thuc nghiem'

read*,tentaptin3

open(333,file=tentaptin3)

do 3333 i=1,M

read(333,'(a4)') datfile(i)

3333 enddo

do 3331 j=1,M

call extract(j,datfile,ia)

do 30 k=1,8195

30 nia(j,k)=ia(j,k)

3331 enddo

open(3333,file='allspecs',status='unknown')

do 300 k=1,8195

300 write(3333,'(25i10)') (nia(j,k),j=1,M)

do 3000 i=1,M

write(ordnum,'(i0)') i

open(33333,file='extru'//ordnum,status='unknown')

do 31 k=1,8195

31 write(33333,'(i25)') nia(i,k)

3000 close(33333)

goto 10

endif

else

print*,'Ban muon xu ly bao nhieu file?'

read*,L3

if(L3.eq.1) then

print*,’Nhap ten file thuc nghiem’

read*,ten4

open(30,file=ten4,form='formatted')

print*,’Nhap ten file xuat ra’

read*,ten5

open(300,file=ten5)

aa=':'

do i=1,14

read(30,'(76a1)')(dummy(j),j=1,76)

enddo

i=0

303 if (.not.eof(30)) then

read(30,'(6a1)') (dummy(j),j=1,6)

if(dummy(6).ne.aa) goto 303

backspace (30)

read(30,330) dumm,(ic(i+l), l=1,8)

i=i+8

max=i

goto 303

else

do i=1,max

ec(i) = 0.0062 + (1.9485-0.0062)/8192*(i-1)

write(300,3003) ec(i), ic(i)

enddo

goto 10

endif

3330 format(a1)

330 format(a6,8i8)

3003 format(f12.6,8i8)

else

print*,'Nhap file chua ten cac file thuc nghiem'

read*,tentaptin33

open(7,file=tentaptin33)

do 7 i=1,L3

read(7,'(a6)') datfile2(i)

7 enddo

do 77 j=1,L3

call extract2(j,datfile2,id)

do 70 k=1,8192

70 nid(j,k)=id(j,k)

77 enddo

open(8,file='allspecs2',status='unknown')

do 700 k=1,8192

700 write(8,'(25i10)') (nid(j,k),j=1,L3)

do 770 i=1,L3

write(ordnum2,'(i0)') i

open(9,file='extruBM'//ordnum2,status='unknown')

do 707 k=1,8192

707 write(9,'(i25)') nid(i,k)

770 close(9)

goto 10

endif

endif

stop UChương trình con

1. Đoạn chương trình xử lý Phổ Trung tâm Hạt nhân TPHCM

subroutine extract(j,datfile,ia)

character*4 datfile(100)

character*1 dummy(76), aa

character*6 dumm

dimension ia(100,10000),e(10000)

open(3,file=datfile(j),form='formatted')

aa=':'

2 do i=1,9

read(3,'(76a1)')(dummy(l),l=1,76)

enddo

i=0

11 if (.not.eof(3)) then

read(3,'(6a1)') (dummy(l),l=1,6)

if(dummy(6).ne.aa) goto 11

backspace (3)

read(3,111) dumm,(ia(j,i+m), m=1,5)

i=i+5

max=i

goto 11

endif

101 format(a1)

111 format(a6,5i13)

301 format(f12.6,5i13)

close(3)

return

end

2. Đoạn chương trình xử lý Phổ của Phòng thí nghiệm Bộ môn Vật lý hạt nhân Đại học Khoa

học tự nhiên TPHCM

subroutine extract2(j,datfile2,id)

character*6 datfile2(100)

character*1 dummy(76), aa

character*6 dumm

dimension id(100,10000)

open(6,file=datfile2(j),form='formatted')

aa=':'

6 do i=1,14

read(6,'(76a1)')(dummy(l),l=1,76)

enddo

i=0

66 if (.not.eof(6)) then

read(6,'(6a1)') (dummy(l),l=1,6)

if(dummy(6).ne.aa) goto 66

backspace (6)

read(6,660) dumm,(id(j,i+m), m=1,8)

i=i+8

max=i

goto 66

endif

606 format(a1)

660 format(a6,8i8)

666 format(f12.6,8i8)

close(6)

return

end

UCửa sổ chương trình

Hình 3.17: Cửa sổ chương trình “Trích xuất dữ liệu từ file thực nghiệm”

Chọn phổ thực nghiệm của Trung tâm Hạt nhân TPHCM hay Phòng thí nghiệm Bộ môn Vật lý

hạt nhân Đại học Khoa học tự nhiên TPHCM bằng cách nhập 1 hoặc 2. Các thao tác sau đó hoàn toàn

tương tự nhau ở cả hai trường hợp. Trong luận văn này, chúng tôi chỉ giới thiệu cách trích xuất dữ liệu

từ file thực nghiệm của Trung tâm Hạt nhân TPHCM.

Hình 3.18: Cửa sổ chương trình “Trích xuất dữ liệu từ file thực nghiệm”- Nhập số file cần xử lý

Nhập số file thực nghiệm muốn xử lý.

Hình 3.19: Cửa sổ chương trình “Trích xuất dữ liệu từ file thực nghiệm”-Nhập tên file thực

nghiệm”

Nếu chỉ xử lý một file thì nhập tên file đó.

Hình 3.20: Cửa sổ chương trình “Trích xuất dữ liệu từ file thực nghiệm”-Nhập tên file xuất ra”

Nhập tên file kết quả xử lý xuất ra (nếu muốn file có định dạng Excel thì chọn phần mở rộng

là .xls).

Trong file thực nghiệm ban đầu, các số liệu được sắp xếp thành nhiều cột, gây bất tiện cho

quá trình tính toán và vẽ phổ gamma. Do đó trong file kết quả thu được các số liệu sẽ được sắp xếp

lại thành hai cột. Cột thứ nhất là kênh năng lượng, cột thứ hai là số đếm tương ứng.

Hình 3.21: Cửa sổ chương trình “Trích xuất dữ liệu từ file thực nghiệm”-Nhập file chứa tên các

file thực nghiệm”

Trong trường hợp người sử dụng cần trích xuất dữ liệu thực nghiệm từ nhiều file, chương

trình sẽ hỏi file chứa tên các file thực nghiệm. Ở ví dụ này, người sử dụng viết tên 3 file thực

nghiệm cần xuất dữ liệu vào một file (xuống dòng để ngăn cách) và đặt tên là thucnghiem. Chương

trình sẽ xuất ra các file kết quả extru1, extru2, extru3 tương ứng.

3.2.4. Modun 4 - Chương trình chuẩn phổ mô phỏng và phổ thực nghiệm

Đối với người sử dụng các chương trình mô phỏng nói chung cũng như MCNP5 nói riêng,

việc so sánh phổ mô phỏng và phổ thực nghiệm là một bước không thể thiếu. Bước này nhằm kiểm

tra tính hiệu lực, độ tin cậy của chương trình mô phỏng trong bài toán về phổ gamma cũng như

kiểm tra độ tin cậy của các thông tin do nhà sản xuất cung cấp. Việc mô phỏng được thực hiện với

cấu hình hệ phổ kế tương tự như thực nghiệm.

Trong phổ thực nghiệm, khi bức xạ đi vào detector thì sẽ ghi nhận được các xung điện thế có

biên độ tỷ lệ với năng lượng của bức xạ tới trong vùng nhạy. Sau khi qua bộ tiền khuếch đại và

khuếch đại chính đến bộ đếm MCA các xung có biên độ khác nhau được sắp xếp vào các kênh khác

nhau. Xung có biên độ cao hơn được xếp vào vị trí kênh cao hơn, trong file kết quả thực nghiệm, ở

mỗi kênh sẽ thể hiện số đếm cho biết có bao nhiêu xung năng lượng có biên độ ứng với kênh này.

Tuy nhiên trong phổ mô phỏng, năng lượng ghi nhận được hiển thị dưới dạng độ cao xung vi phân.

Vì vậy, để có thể so sánh phổ mô phỏng và phổ thực nghiệm, cần một thao tác gọi là “chuẩn phổ”.

Trong quá trình chuẩn phổ này, chúng tôi tiến hành chuyển phân bố độ cao xung vi phân trong file

mô phỏng thành dạng số đếm tương ứng với file thực nghiệm bằng cách chọn một đỉnh năng lượng

làm chuẩn.

Đỉnh năng lượng được chọn để chuẩn phổ thường là đỉnh năng lượng có giá trị lớn nhất trong

phổ gamma. Chọn lựa này nhằm hạn chế tối đa ảnh hưởng bởi các nền Compton của các quang đỉnh

khác.

USơ đồ khối chương trình

Bắt đầu

Nhập tên file trích xuất dữ liệu từ output

Nhập tên file trích xuất dữ liệu từ thực nghiệm

Nhập đỉnh năng lượng chuẩn phổ

Nhập tên file kết quả xuất ra

Mở file output đã sắp xếp, xác định vị trí, lấy ra các giá trị độ cao xung năng lượng

Mở file thực nghiệm, lấy ra các giá trị số đếm

Xác định vị trí ứng với năng lượng chuẩn phổ đã được nhập vào. Rút ra xung năng lượng vi phân tương ứng

Tính tỷ số giữa số đếm và giá trị xung năng lượng vi phân

Chuyển các giá trị độ cao xung sang số đếm bằng cách

nhân mỗi giá trị xung năng lượng vi phân với tỷ số trên

Hiển thị file kết quả

Kết thúc

Hình 3.22: Sơ đồ khối chương trình “Chuẩn phổ”

UĐoạn chương trình

Program Chuanpho

use portlib

use msflib

use dflib

character*10 outfile,thucnghiem,kqchuanpho

character*1 dum(200)

real*4 zjcount(10000),de4(10000),ratio,zde(10000),

+e4(10000),min4,hieu4(10000),dinhnlcp

integer jcount(10000),dinh

print*,'****************************************************'

print*,'*CHUONG TRINH CHUAN PHO MO PHONG VA PHO THUC NGHIEM*'

print*,'****************************************************'

print*,''

print*,'Nhap ten file trich xuat du lieu tu output '

read*,outfile

print*,'Nhap ten file trich xuat du lieu tu thuc nghiem'

read*,thucnghiem

print*,'Nhap dinh nang luong chuan pho(keV)'

read*,dinhnlcp

print*,'Nhap ten file ket qua xuat ra'

read*,kqchuanpho

open(4,file=outfile)

do 46 k=1,9

46 read(4,'(40a1)')(dum(l),l=1,40)

do 4 k=1,8194

4 read(4,'(16a1,en12.3,65a1)')(dum(l),l=1,16),de4(k),

+(dum(l),l=1,65)

open(40,file=thucnghiem)

do 40 k=1,8194

40 read(40,43)(dummy(l),l=1,20),jcount(k),

+(dummy(l),l=1,65)

43 format(20a1,i5,65a1)

do 415 k=1,8194

415 hieu4(k)=abs(e4(k)*1000-dinhnlcp)

min4=hieu4(1)

dinh=1

do 405 k=2,8194

if (min4.gt.hieu4(k)) then

min4=hieu4(k)

dinh=k

405 endif

ratio=jcount(dinh)/de4(dinh)

do 41 k=1,8194

41 de4(k)=de4(k)*ratio

open(41,file=kqchuanpho,status='unknown')

write(41,48)

48 format('Dinh nang luong chuan pho (keV):')

write(41,'(f12.4)') e4(dinh)*1000

write(41,47)

47 format('Kenh nang luong(MeV) So dem thuc nghiem So dem + mo

phong')

do 42 k=1,8194

42 write(41,'(f12.7,10x,i5,10x,f12.0)')e4(k),jcount(k),de4(k)

close (41)

goto 10

stop

UCửa sổ chương trình

Hình 3.23: Cửa sổ chương trình “Chuẩn phổ” - Nhập tên file trích xuất dữ liệu từ output

Nhập tên file output đã trích xuất và sắp xếp dữ liệu, file này là kết quả của modun 2: “Trích

xuất dữ liệu từ file output”.

Hình 3.24: Cửa sổ chương trình “Chuẩn phổ” - Nhập tên file trích xuất dữ liệu từ thực nghiệm

Nhập tên file thực nghiệm đã trích xuất và sắp xếp dữ liệu, file này là kết quả của modun 3:

“Trích xuất dữ liệu từ file thực nghiệm”

Hình 3.25: Cửa sổ chương trình “Chuẩn phổ”- Nhập đỉnh năng lượng chuẩn phổ

137

PCs, đỉnh năng lượng ở mức 661,66 keV.

Chọn đỉnh năng lượng muốn chuẩn phổ, lấy đơn vị là keV. Đó thường là đỉnh năng lượng có

giá trị lớn nhất của phổ. Đối với P

Hình 3.26: Cửa sổ chương trình “Chuẩn phổ”- Nhập tên file kết quả xuất ra

Nhập tên file kết quả chuẩn phổ xuất ra. Trong file này, ứng với mỗi kênh năng lượng sẽ là

số đếm tương ứng của phổ thực nghiệm và phổ mô phỏng, từ đó, người sử dụng dễ dàng xác định

độ tin cậy của chương trình mô phỏng dựa trên so sánh hai số đếm này. Kết quả này cũng được sử

dụng để vẽ đồ thị biểu diễn phổ gamma mô phỏng và phổ gamma thực nghiệm. Đây là hình ảnh so

sánh rất trực quan và sẽ được nghiên cứu kỹ hơn trong modun 7 của chương trình.

3.2.5. Modun 5 - Chương trình tính diện tích đỉnh phổ và hiệu suất

Một trong những vấn đề thường được nghiên cứu là tính diện tích đỉnh phổ năng lượng và

hiệu suất của quá trình thực nghiệm hoặc quá trình mô phỏng. Chúng tôi thiết kế modun này nhằm

tính toán hiệu suất mô phỏng dựa trên những số liệu trích xuất được từ file output. Trong modun

này, chúng tôi chia ra hai trường hợp, tính hiệu suất đối với phổ gamma đơn năng (chỉ có một đỉnh

năng lượng) và tính hiệu suất đối với phổ gamma đa năng (chỉ xét trường hợp có hai đỉnh năng

lượng trong phổ).

Đối với một hệ phổ kế gamma cụ thể, bố trí hình học đo xác định và tại vạch năng lượng

quan tâm thì hiệu suất detector có giá trị xác định. Do đó, hiệu suất là một trong những thông số

quan trọng dùng để nghiên cứu các đặc trưng của detector, nguồn phóng xạ và hình học đo.

Trong đo đạc thực nghiệm trên hệ phổ kế gamma, khái niệm “hiệu suất” được hiểu là hiệu

suất đỉnh năng lượng toàn phổ (full energy peak effciency) được định nghĩa là tỷ số giữa tốc độ đếm

đỉnh ở năng lượng E (số đếm chia cho thời gian đo) và tốc độ phát gamma từ nguồn cũng ở năng

lượng E tương ứng.

Trong tính toán MCNP5, hiệu suất của detector được xác định bằng tỷ số giữa số photon phát

ra từ nguồn theo mọi hướng và số photon đóng góp vào quang đỉnh của phổ gamma mô phỏng.

Thực tế khi tính hiệu suất đỉnh năng lượng toàn phần của phổ gamma mô phỏng còn phải xét đến

một số yếu tố như tỷ số giữa số hạt để lại năng lượng khác 0 trong detector và tổng số hạt khảo sát,

suất phát tương ứng với mỗi đỉnh năng lượng trong trường hợp tính hiệu suất của phổ gamma đa

năng,…

USơ đồ khối chương trình

Nhập lựa chọn Bắ 2 Phổ đơn phổ gamma đơn t 1

Nhập giá trị đỉnh

Nhập suất phát tương ứng với đỉnh năng lượng 2 năng lượng (keV)

1

2 Nhập bề rộng đỉnh 1 2

Nhập tên file trích xuất

dữ liệu từ output

1

2 Nhập tên file kết quả xuất ra

1 2

Mở file trích xuất dữ liệu từ output, lấy ra

các giá trị độ cao xung năng lượng vi phân

1 2

Mở file trích xuất dữ liệu từ output, lấy ra các giá trị

tổng số hạt khảo sát, tổng số hạt để lại năng lượng

ổ

2 1 Kết nối chương

2 trình con tính

diện tích đỉnh Tính hiệu suất bằng tỷ số giữa diện tích đỉnh năng lượng trên tổng độ cao xung năng lượng vi phân, nhân với tỷ số giữa số hạt bỏ năng lượng khác 0 trong detector trên tích của tổng số hạt khảo sát và suất phát

1

Tính hiệu suất bằng tỷ số giữa diện tích đỉnh năng lượng trên tổng 2 độ cao xung năng lượng vi phân, nhân với tỷ số giữa số hạt bỏ

ổ ố

Hiển thị file kết quả Kết thúc

Hình 3.27: Sơ đồ khối chương trình “Tính diện tích đỉnh phổ và hiệu suất”

1. Đúng; 2. Sai

USơ đồ khối của chương trình con

Bắt đầu

Xác định bề rộng vùng đỉnh năng lượng

Xác định vùng xung quanh đỉnh năng lượng để tính năng lượng phông

Xác định giá trị phông năng lượng xung quanh đỉnh

Tính diện tích đỉnh năng lượng (bao gồm phông)

Tính diện tích đỉnh năng lượng đã trừ phông

Kết thúc

UĐoạn chương trình

Hình 3.28: Sơ đồ khối chương trình con “Tính diện tích đỉnh phổ và hiệu suất”

Program Tinh dien tich dinh pho va hieu suat

use portlib

use msflib

use dflib

integer phpos,phwid,l,r,lside,rside,NL,nht5,nps5

real*8 de5(8194),lsideba,rsideba,tpa,ba,

+ eff,pa,pc,nat5,e5(8194),dinhnl,min,hieu(8194),sp

character*10 outfile5,kqhieusuat

print*,'*************************************************'

print*,'*CHUONG TRINH TINH DIEN TICH DINH PHO VA HIEU SUAT*'

print*,'*************************************************'

print*,''

print*,'Pho nang luong:'

print*,'1- Don nang (mot dinh nang luong)'

print*,'2- Da nang (co hai dinh nang luong)'

read*,NL

select case (NL)

case(1)

print*,'nhap gia tri dinh nang luong (keV)'

read*,dinhnl

print*,'nhap be rong dinh'

read*,phwid

print*,'nhap ten file trich xuat du lieu tu output'

read*,outfile5

print*,'Nhap ten file ket qua xuat ra:'

read*,kqhieusuat

open(5,file=outfile5)

do 51 k=1,3

51 read(5,'(50a1)') (dum(l),l=1,50)

read(5,'(f12.4)') nat5

read(5,'(50a1)') (dum(l),l=1,50)

read(5,'(i12)') nps5

read(5,'(50a1)') (dum(l),l=1,50)

read(5,'(i12)') nht5

read(5,'(50a1)') (dum(l),l=1,50)

do 5000 k=1,8193

5000 read(5,'(f12.4,5a1,en14.3,65a1)')e5(k),

+(dum(l),l=1,5),de5(k),(dum(l),l=1,65)

do 515 k=1,8193

515 hieu(k)=abs(e5(k)*1000-dinhnl)

min=hieu(1)

phpos=1

do 505 k=2,8193

if (min.gt.hieu(k)) then

min=hieu(k)

phpos=k

505 endif

call compupa(phpos,phwid,de5,pa)

eff=pa/(nat5-de5(1)-de5(2))*nht5/nps5

pc=pa/(nat5-de5(1)-de5(2))*nht5

case(2)

print*,'nhap gia tri dinh nang luong (keV)'

read*,dinhnl

print*,'nhap suat phat tuong ung voi dinh nang luong'

read*,sp

print*,'nhap be rong dinh'

read*,phwid

print*,'nhap file output da sap xep'

read*,outfile5

print*,'Nhap ten file ket qua xuat ra:'

read*,kqhieusuat

open(5,file=outfile5)

do 512 k=1,3

512 read(5,'(50a1)') (dum(l),l=1,50)

read(5,'(f12.4)') nat5

read(5,'(50a1)') (dum(l),l=1,50)

read(5,'(i12)') nps5

read(5,'(50a1)') (dum(l),l=1,50)

read(5,'(i12)') nht5

read(5,'(50a1)') (dum(l),l=1,50)

do 5002 k=1,8193

5002 read(5,'(f12.4,5a1,en14.3,65a1)')e5(k),

+(dum(l),l=1,5),de5(k),(dum(l),l=1,65)

do 5152 k=1,8193

5152 hieu(k)=abs(e5(k)*1000-dinhnl)

min=hieu(1)

phpos=1

do 5052 k=2,8193

if (min.gt.hieu(k)) then

min=hieu(k)

phpos=k

5052 endif

call compupa(phpos,phwid,de5,pa)

eff=pa/(nat5-de5(1)-de5(2))*nht5/(nps5*sp)

pc=pa/(nat5-de5(1)-de5(2))*nht5

end select

open(50,file=kqhieusuat)

write(50,555)

555 format(' dinh NL(keV) dien tich dinh hieu suat')

write(50,'(f12.3,5x,f12.7,4x,f12.7)')1000*e5(phpos),pa,eff

close(5)

close(50)

goto 10

stop

UChương trình con

subroutine compupa(phpos,phwid,denew,pa)

integer phpos,phwid,l,r,lside,rside

real*8 denew(8194),lsideba,rsideba,tpa,ba,pa

l=phpos-phwid

r=phpos+phwid

lside=phpos-(2*phwid-9)

rside=phpos+(2*phwid-9)

lsideba=0

do 10 i=lside,l-1

10 lsideba=lsideba+denew(i)

rsideba=0

do 20 i=r+1,rside

20 rsideba=rsideba+denew(i)

tpa=0

do 30 i=l,r

30 tpa=tpa+denew(i)

ba=(lsideba+rsideba)*phwid/(phwid-9)

pa=tpa-ba

return

end UCửa sổ chương trình

Hình 3.29: Cửa sổ chương trình “Tính diện tích đỉnh phổ và hiệu suất”- Chọn phổ năng lượng

Nhập 1 nếu phổ gamma là phổ đơn năng (chỉ có một đỉnh năng lượng).

Hình 3.30: Cửa sổ chương trình “Tính diện tích đỉnh phổ và hiệu suất”- Nhập giá trị đỉnh năng

lượng

Nhập giá trị đỉnh năng lượng cần tính diện tích và hiệu suất, đơn vị là keV.

Hình 3.31: Cửa sổ chương trình “Tính diện tích đỉnh phổ và hiệu suất”- Nhập bề rộng đỉnh năng

lượng

Nhập bề rộng đỉnh năng lượng cần tính diện tích, bề rộng đỉnh thường chọn là 12 (khoảng 12

kênh năng lượng), bề rộng này là vừa đủ để tính diện tích nền đỉnh quan tâm mà không bị ảnh

hưởng bởi những đỉnh năng lượng khác liền kề.

Hình 3.32: Cửa sổ chương trình “Tính diện tích đỉnh phổ và hiệu suất”- Nhập tên file trích xuất dữ

liệu từ output

Nhập tên file output đã được trích xuất và sắp xếp dữ liệu (file kết quả thu được từ modun 2).

Hình 3.33: Cửa sổ chương trình “Tính diện tích đỉnh phổ và hiệu suất”- Nhập tên file kết quả xuất

ra

Nhập tên file kết quả xuất ra.

Trong trường hợp phổ gamma có hai đỉnh năng lượng, chương trình sẽ yêu cầu nhập thêm

suất phát tương ứng của mỗi đỉnh. Số liệu này lấy từ thư viện số liệu vật lý hạt nhân.

Hình 3.34: Cửa sổ chương trình “Tính diện tích đỉnh phổ và hiệu suất”- Nhập suất phát tương ứng

với đỉnh năng lượng

Khi nhập suất phát, nên nhập con số với độ chính xác cao như trong hình trên vì trong lúc

tính toán hiệu suất từng đỉnh năng lượng sẽ phải tính tích của tổng số hạt đầu vào để khảo sát và

suất phát, tổng số hạt đầu vào để khảo sát thường có giá trị rất lớn (vài trăm triệu đến vài trăm tỷ),

do đó một sự thay đổi nhỏ trong suất phát cũng làm thay đổi đáng kể hiệu suất tính được.

2.3.6. Modun 6 - Chương trình vẽ hệ phổ kế

Modun này có chức năng kết nối chương trình Fortran với phần mềm vẽ hệ phổ kế của

MCNP5. Khi lựa chọn modun này, cửa sổ chương trình Vised của MCNP5 sẽ xuất hiện, việc vẽ phổ

UĐoạn chương trình

sẽ được tiến hành trên cửa sổ này.

Program Vehephoke

logical result6

result6 = systemqq('D:\vised')

UCửa sổ chương trình

goto 10

Hình 3.35: Cửa sổ chương trình “Vẽ hệ phổ kế”

Trong cửa sổ này, người sử dụng mở file input cần vẽ hệ phổ kế rồi chọn Update, sau đó sẽ

được kết quả như hình 3.36.

Hình 3.36: Cửa sổ chương trình “Vẽ hệ phổ kế”- hình ảnh mặt cắt dọc hệ phổ kế

Có nhiều chức năng cho người sử dụng lực chọn để chỉnh sửa hình ảnh theo mong muốn.

2.3.7. Modun 7 - Hình ảnh phổ gamma

UĐoạn chương trình

Modun này sẽ kết nối chương trình MCNP5 với một số phổ gamma dựng sẵn.

Program Hinhanhphogamma

logical result3

result3 = systemqq('D:\Phogamma.xls')

goto 10

Khi lựa chọn modun này, người sử dụng sẽ nhận được một file định dạng Excel, file này

60

137

54

232

238

241 (P

PAm, P

PCo, P

PCs, P

PMn, P

PTh, P

PU). Đây là các nguồn điểm đặt dọc detector và cách detector

chứa dữ liệu phổ mô phỏng của một số nguồn phóng xạ thường gặp

15cm. Các dữ liệu này là kết quả thu được sau quá trình tiến hành mô phỏng (trên hệ phổ kế của

Trung tâm hạt nhân TPHCM) và xử lý phổ gamma (lập input, chạy MCNP5, xử lý output, tính toán

60

PCo

diện tích đỉnh năng lượng, hiệu suất, vẽ phổ gamma…).

2000

Hình sau minh họa phổ gamma mô phỏng của P

1332 keV

1600

1200

1173 keV

m ế đ ố S

800

BS

400

0

0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

DE 511 keV SE

Năng lượng (MeV)

60

PCo

60

Hình 3.37 : Phổ gamma mô phỏng của P

PCo là nguồn được sử dụng nhiều trong thực tế, đồng thời năng lượng gamma của

60

PCo có thể gây ra hầu hết các hiệu ứng đại diện cho tương tác của bức xạ gamma với

Nguồn P

nguồn P

detector. Nguồn được đo ở khoảng cách 15 cm so với bề mặt detector trong thời gian 14400 giây, số

lịch sử hạt là 900000000 nhằm loại bỏ những thăng giáng thống kê trong phổ gamma.

Trong phổ gamma trên, trục hoành biểu thị năng lượng, đơn vị MeV, trục tung là số đếm do

detector ghi nhận (dưới dạng phân bố độ cao xung). Trên nền Compton liên tục sẽ có những đỉnh

năng lượng sau:

Đỉnh năng lượng toàn phần ở 1173 keV và 1332 keV được hình thành bởi các tương tác mà ở

đó tia gamma mất toàn bộ năng lượng của nó bằng hiệu ứng quang điện hay một chuỗi các tán xạ

Compton và kết thúc cũng bằng một hiệu ứng quang điện.

Đỉnh nhỏ ở vị trí 310 keV là đỉnh thoát đôi (DE) còn đỉnh nhỏ ở vị trí 830 keV là đỉnh thoát

đơn (SE). Ở giữa hai đỉnh này là đỉnh bức xạ hủy (511 keV). Các đỉnh này hình thành do quá trình

tạo cặp.

CHƯƠNG 4: KẾT LUẬN VÀ KIẾN NGHỊ

KẾT LUẬN CHUNG

Từ những mục tiêu đặt ra ban đầu, chúng tôi đã đạt được những kết quả cụ thể sau đây:

- Đã kết nối được chương trình MCNP5 và chương trình Fortran để có thể chạy tự động các

file input được liệt kê, đây là một thành công quan trọng vì nó là bước đầu tiên để có thể thực hiện

các bước tiếp theo của quá trình dùng ngôn ngữ lập trình Fortran tăng tính hiệu quả trong việc sử

dụng chương trình MCNP5. Đối với mục tiêu này, chúng tôi đã đạt được yêu cầu đặt ra ban đầu,

chương trình có thể chạy tự động một file input hoặc nhiều file input như ý muốn mà không gặp khó

khăn gì.

- Đã xử lý được thông tin từ file output một cách tự động, trích xuất được những dữ liệu cần

thiết một cách đầy đủ và chính xác. Trong phần này, chúng tôi cũng có thể trích xuất dữ liệu từ một

file hoặc nhiều file theo ý muốn. Đối với dữ liệu có được từ các chương trình thực nghiệm, chúng

tôi cũng đã dùng Fortran để xử lý nhanh chóng các file kết quả thực nghiệm từ Trung tâm Hạt nhân

TPHCM và từ Phòng thí nghiệm bộ môn Vật lý hạt nhân Đại học Khoa học tự nhiên TPHCM.

- Đã thiết lập được chương trình chuẩn phổ thực nghiệm và phổ lý thuyết, từ đó có thể so sánh

độ tin cậy của kết quả có được từ chương trình mô phỏng MCNP5, đó cũng là những số liệu cơ sở

để tiến hành vẽ và so sánh phổ gamma mô phỏng và phổ gamma thực nghiệm.

- Xác định được diện tích đỉnh phổ gamma và hiệu suất. Đây là một trong những ứng dụng

tiếp theo của chương trình “AUT-MCNP5”. Diện tích đỉnh phổ và hiệu suất sẽ được tính dựa vào

năng lượng đỉnh, bề rộng đỉnh, suất phát được nhập vào. Chương trình tính được hiệu suất trong cả

hai trường hợp phổ gamma đơn năng và phổ gamma đa năng với hai đỉnh năng lượng.

- Kết nối nhanh với ứng dụng vẽ hình của MCNP5 bằng modun vẽ hình hệ phổ kế. Đây là

modun biểu diễn hình ảnh trực quan, giúp người sử dụng chương trình có thể tiện lợi trong việc mô

hình hóa hệ phổ kế bao gồm thiết lập, sửa chữa và cập nhật input MCNP5. Chương trình Fortran chỉ

có nhiệm vụ gọi ra cửa sổ vẽ hình này, còn những thao tác tiếp theo được tiến hành trực tiếp trên

cửa sổ chương trình Vised của MCNP5.

- Đã xây dựng modun truy xuất một số phổ gamma điển hình. Phần hình ảnh phổ gamma là

kết quả sau khi đã phân tích và trích xuất dữ liệu từ output.

KIẾN NGHỊ HƯỚNG NGHIÊN CỨU TIẾP THEO

Chương trình AUT-MCNP5 đã được phát triển và đạt được một số kết quả nhất định giúp

cho người sử dụng MCNP5 trong lĩnh vực mô phỏng phổ gamma giảm bớt các thao tác xử lý, tiết

kiệm thời gian, nâng cao độ tin cậy của các kết quả tính toán. Tuy nhiên, thời gian thực hiện của

luận văn có hạn nên chương trình AUT-MCNP5 vẫn còn một số hạn chế:

- Các modun hiện có của chương trình MCNP5 chỉ mới giải quyết một số vấn đề cơ bản của

quá trình mô phỏng phổ gamma một cách tổng quát nhất. Còn nhiều vấn đề khác như: “thiết lập

input tự động”, “tính toán và phân tích biến đổi hiệu suất theo những yếu tố khác nhau”…chúng tôi

chưa có điều kiện xử lý.

- Việc sử dụng các modun của chương trình AUT-MCNP5 vẫn còn phải tiến hành riêng lẻ

từng modun.

- Cửa sổ chương trình AUT-MCNP5 là cửa sổ DOS, giao diện chưa được bắt mắt và thuận

tiện.

Do đó, để hoàn thiện chương trình MCNP5 nên tiếp tục bổ sung:

- Một số modun có thể được nghiên cứu sâu hơn nhằm xử lý được vấn đề một cách tổng quát

hơn. Bên cạnh đó, người nghiên cứu tiếp theo có thể tiếp tục xây dựng những modun khác, xử lý

những vấn đề khác.

- Một hướng khác được đặt ra là hoàn thiện chương trình để có thể tạo nên một chương trình

xử lý liền mạch hoàn chỉnh, với những yêu cầu cụ thể cùng với file input ban đầu, chương trình sẽ

tiến hành tự động trong tất cả các bước tiếp theo và xuất ra những kết quả như mong muốn của

người sử dụng.

- Cũng có thể thiết kế giao diện đẹp mắt, dễ sử dụng hơn và nâng cấp thành một phần mềm

ứng dụng nhiều chức năng dành cho người cần sử dụng MCNP5 trong nghiên cứu.

TÀI LIỆU THAM KHẢO

UTiếng Việt

[1] Ngô Quang Huy (2006). Cơ sở vật lý hạt nhân. Nhà xuất bản Khoa học và kỹ thuật.

[2] Ngô Quang Huy, Đỗ Quang Bình, Võ Xuân Ân (2006). Mô phỏng các phổ gamma phức tạp đo

trên hệ phổ kế gamma dùng detector HPGe bằng chương trình MCNP. Tạp chí phát triển KH &

CN, ĐHQG TPHCM, Tập 9, Số 9, Trang 63-70.

[3] Ngô Quang Huy, Trần Văn Luyến và Nguyễn Văn Mai (1999). Khảo sát nền phông phóng xạ

đối với một số đặc trưng môi trường tại TP Hồ Chí Minh. Báo cáo kết quả nghiên cứu đề tài cấp

Bộ các năm 1996 – 1998, Trung tâm Hạt nhân TPHCM, Viện NLNT Việt Nam.

[4] Trần Văn Luyến (2005). Nghiên cứu nền phông phóng xạ vùng Nam bộ Việt Nam. Luận án tiến

sỹ, Trường Đại học Khoa học Tự nhiên. Đại học Quốc gia TPHCM.

[5] Trương Thị Hồng Loan (2010). Áp dụng phương pháp mô phỏng Monte Carlo để nâng cao chất

lượng hệ phổ kế gamma sử dụng detector bán dẫn HPGe. Luận án tiến sĩ, Trường Đại học Khoa

học Tự nhiên, Đại học Quốc gia TPHCM.

[6] Võ Văn Hoàng (2004). Mô phỏng trong Vật lý. Nhà xuất bản Đại học Quốc gia TPHCM.

[7] Võ Văn Hoàng (2007). Ngôn ngữ lập trình Fortran. Nhà xuất bản Giáo dục.

[8] Võ Xuân Ân (2008). Nghiên cứu hiệu suất ghi nhận của detector bán dẫn siêu tinh khiết (HPGe)

trong phổ kế gamma bằng phương pháp Monte Carlo và thuật toán di truyền. Luận án tiến sĩ,

Trường Đại học Khoa học Tự nhiên, Đại học Quốc gia TPHCM.

[9] Võ Xuân Ân, Mai Văn Nhơn (2002). Phương pháp Monte Carlo trong mô phỏng tương tác của

bức xạ gamma với vật chất. Hội nghị Khoa học Trường Đại học Khoa học tự nhiên lần thứ III,

UTiếng Anh

ĐHQG TPHCM, Tóm tắt, Trang 99. [8]

[10] Canberra Industries, Inc (1995). Ultra Low Background Detector Systems. Canberra

P Edition.

th [11] Canberra Industries, Inc (2000). Germanium Detectors - User's Manual, 12P

Industries, Inc., Connecticut.

Canberra Industries, Inc., Connecticut.

[12] DIGITAL Visual Fortran Language Reference (1998). DIGITAL Visual Fortran Version 6.0.

Standard and Professional Editions.

[13] G. Haase, D.Tait, A. Wiechen (1995). Determination of full energy peak efficiency for

cylindrical volume sources by the use of a point source standard in gamma spectrometry.

Nucl. Instrum. and Methods in Phys. Res. A361, 240-244.

[14] J. F. Briesmeister (2000). MCNP – A General Purpose Monte Carlo N – Particle Transport

Code. Version 4C2, Los Alamos, LA.

[15] Knoll G.F. (1999). Radiation Detection and Measurement, Third Edition. John Wiley & Sons,

Inc., New York.

[16] M.A. Ludington, R.G. Helmer (2000). High accuracy measurememt and Monte Carlo

calculations of the relative efficiency curve of an HPGe detector from 433 to 2754 keV. Nucl.

Instrum. and Methods in Phys. Res. A446, 506-521.

[17] R.D.Evans (1955). The Atomic Nuclear.

[18] X-5 Monte Carlo Team (2003). MCNP - A General Purpose Monte Carlo N-Particle Transport

Code, Version 5, Volume I: Overview and Theory. Los Alamos National Laboratory, LA-UR-

03-1987.

PHỤ LỤC

UPhụ lục 1U: Ảnh chụp hệ phổ kế gamma phông thấp của Trung tâm Hạt nhân TPHCM

UPhụ lục U2: Một input điển hình của chương trình MCNP5

1

2

3

4

5

6

7

8

9

10

11

Problem - HPGe coaxial detector efficiencies and pulse height distribution c Cell cards 1 6 -8.94 -1 -23 21 imp:p,e=1 2 1 -5.35 (-55 -64 22)#1 imp:p=1 imp:e=0 $ cell detector 3 3 -0.00129 (1 -2 -22 21)#1 imp:p,e=1 4 2 -2.6989 (2 -3 -24 21):(-3 -21 20) imp:p,e=1 5 3 -0.00129 (-4 -25 24):(3 -4 -24 20):(-4 -20 16) imp:p,e=1 6 2 -2.6989 (-5 -26 25):(4 -5 -25 16):(-5 -16 15) imp:p,e=1 11 3 -0.00129((-13 -31 30):(-11 -30 26):(5 -9 -26 19)) - - - - - - - - - -#46 imp:p,e=1 $ - cs40e

12 6 -8.94 9 -10 -26 19 imp:p,e=1 13 7 -0.88 (11 -12 -30 26):(10 -12 -26 19) imp:p,e=1 - -

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

- - - - - - - - - - - - -

14 8 -7.28 12 -13 -30 19 imp:p,e=1 15 9 -11.34 (13 -14 -31 17):(5 -14 -17 16) imp:p,e=1 16 3 -0.00129 (-13 -32 31):(13 -14 -34 31) imp:p,e=1 17 8 -7.28 -13 -33 32 imp:p,e=1 18 6 -8.94 -13 -34 33 imp:p,e=1 19 10 -7.86 13 -14 -35 34 imp:p,e=1 20 9 -11.34 (-14 -36 35):(-13 -35 34) imp:p,e=1 21 10 -7.86 5 -14 -16 15 imp:p,e=1 22 8 -7.28 5 -13 -18 17 imp:p,e=1 23 6 -8.94 5 -13 -19 18 imp:p,e=1 24 0 14:36:-15 imp:p,e=0 34 15 -2.31 (1 -54 -23 22):(-54 -56 23) imp:p,e=1 $ dien cuc cay ion B 35 17 -5.05 (-55 -24 64):(55 -2 -24 22) imp:p,e=1 $ dien cuc khuech tan Li - th22a25 36 18 -1.435 -3 -81 24 imp:p,e=1 $ cua so IR mylar - cs40e 45 12 -8.92 -76 -79 78 imp:p,e=1 $ nguon Co - co40e 46 13 -1.15 (-77 -80 78)#45 imp:p,e=1 $ holder epoxy - cs40e 47 20 -1.11 -3 -57 81 imp:p,e=1 $ cua so IR kapton - cs40e c Surface cards 1 cz 0.35 2 cz 2.7 $ cs36a 3 cz 2.93 4 cz 3.66 5 cz 3.81 9 cz 7.35 10 cz 7.95 11 cz 9.45 12 cz 14.2 13 cz 15.0 14 cz 25.0 15 pz 0.0 16 pz 1.6 17 pz 10.0 18 pz 10.8 19 pz 11.6 20 pz 19.715 21 pz 19.815 22 pz 20.815 23 pz 22.515 24 pz 24.015 - - - - - - - - - - - - - - - - - - - - - - - - - -

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

25 pz 24.55 26 pz 24.7 30 pz 35.8 31 pz 43.6 32 pz 44.1 33 pz 44.5 34 pz 44.6 35 pz 46.2 36 pz 54.2 54 cz 0.3503 $ dien cuc khuech tan Li - i42 55 cz 2.665 $ dien cuc khuech tan Li - cs36a 56 pz 22.5153 $ dien cuc loi B - i42 57 pz 24.026 $ cua so IR - cs36a 64 pz 23.98 $ dien cuc cua so Li - cs36a 76 cz 0.05 $ hinh hoc nguon Co (2) - co40e 77 cz 1.27 $ hinh hoc nguon Co (2) - co40e 78 pz 39.7 $ hinh hoc nguon Co (2) - co40f 79 pz 39.8 $ hinh hoc nguon Co (2) - co40f 80 pz 40.34 $ hinh hoc nguon Co (2) - co40f 81 pz 24.016 $ split mylar and kapton - cs40e c Data cards mode p m1 32000 -1.0 $ Ge m2 13000 -1.0 $ Al m3 7000 -0.755 8000 -0.232 18000 -0.013 $ Atmosphere m6 29000 -1.0 $ Cu m7 1000 -0.1549 6000 -0.8451 $ Paraffin C9H20 m8 50000 -1.0 $ Sn m9 82000 -1.0 $ Pb m10 26000 -1.0 $ Fe m12 27000 -1.0 $ Co m13 1000 -0.06 6000 -0.721 8000 -0.219 $ Epoxy - cs40e m15 5000 -1.0 $ B m17 32000 -0.9 3000 -0.1 $ Ge Li - th22a24 m18 1000 -0.053 6000 -0.526 8000 -0.421 $ Mylar C10H12O6 - cs40e m20 1000 -0.028 6000 -0.720 7000 -0.077 8000 -0.175 $ Kapton C22H10N2O4 - cs40e sdef cel=45 pos=0 0 0 axs=0 0 1 ext=d1 rad=d2 erg=d3 par=2 wgt=10 ft8 geb 0.00071 0.00075 0.46493 si1 h 39.7 39.8 - - -

90

91

92

93

94

95

96

97

98

99

100

101

102

sp1 d -21 0 si2 h 0.0 0.05 sp2 -21 1 si3 l 1.17324 1.33250 sp3 d 0.999 0.999824 $ co40e f8:p 2 e8 0 .0001 .005471 8190i 1.942341 $ cs40e phys:p $ produce bremsstrahlung radiations - cs40e phys:e cut:p 2j 0 0 $ because of taking a tally of pulse height distributions cut:e nps 900000000 ctme 240 - - - - - - - - - - - - -