Panin_M_P_Modelirovanie_perenosa_izluchenia
.pdf• |
∞ |
∞ |
r |
′ r′ r r r′ |
r * r′ |
′ r′ |
||
|
|
|||||||
J рас = ∫ ∫ ∫ ∫∫ ∫C(E, Ω → E , Ω | r )T (r → r | E,Ω)P (r , E , Ω )× |
||||||||
|
∞ 0 4π ∞ 0 4π |
r |
r |
r |
r |
|
|
|
|
|
|
|
|
||||
|
|
|
×ψ(r |
, E, Ω) dr |
′dE′dΩ′ dr dE dΩ. |
|
||
С |
учетом |
конкретной (сингулярной) |
функции |
детектора |
||||
* |
r′ ′ r′ |
r′ |
r |
r′ |
|
и дельта-функции, содержащейся |
||
P |
(r , E ,Ω ) =δ (r |
−rd )/ |
Σ(r , E) |
|
в транспортном ядре, данный интеграл упрощается:
• |
∞ |
∞ |
|
|
|
J рас = ∫ ∫ ∫ ∫ |
||
|
∞ 0 4π |
0 |
r |
′ |
|
r |
r |
r |
e |
−τ (rr→rrd , E′) |
|
|||
|
rd −r |
|
|
|
|
|
|||||
C(E, Ω → E |
, |
r |
r |
| r ) |
|
|
|
|
|
× |
|
|
|
r |
r |
2 |
|||||||
|
|
|
rd −r |
|
|
|
r |
−r |
|
(11.13) |
|
r |
r |
|
′ |
r |
|
|
d |
|
|
|
|
|
|
|
|
|
|
||||||
×ψ (r, E,Ω) dE dr dE dΩ. |
|
|
|
Для многих видов столкновения между энергией после рассеяния и изменением направления движения частицы существует жесткая связь. Это позволяет провести аналитическое интегрирование и по E′.
Запишем явный вид выражения (11.13) в случае комптоновского рассеяния, дифференциальное макроскопическое
сечение которого есть (см. формулу (10.11)): |
|
|
|||||
r |
′ |
′ |
r |
r′ |
r |
′ |
− Es ). (11.14) |
ΣKN (E, Ω → E , Ω |
|
| r ) = ΣKN (Ω → Ω |
| r, E) ×δ (E |
|
|||
|
|
|
|
|
|
|
′ |
Для компактности обозначим косинус угла рассеяния μs = Ω Ω , а |
|||||||
энергию, выраженную в относительных единицах, α = E / ε . Тогда |
угловую часть дифференциального сечения можно записать в виде
r |
r |
|
r |
r |
2 |
|
|
1 |
|
|
|
2 |
|
|
|
|
||
|
|
′ |
r0 |
|
|
|
|
|
|
|
|
|
|
|
||||
|
|
| r, E,) = ne (r ) |
|
|
|
|
|
|
|
|
× |
|
|
|
||||
ΣKN (Ω → Ω |
2 |
|
|
|
|
|
|
|
|
|
|
|||||||
|
|
|
|
|
1 |
+α(1− μs ) |
|
|
|
|
||||||||
|
|
|
|
|
|
|
1 |
|
|
|
|
|
|
|
|
2 |
|
|
|
|
|
|
|
|
|
|
|
|
+1 |
+α(1 |
− μs ) −(1− μs ) |
|
|
, |
|||
|
|
|
|
|
|
|
|
|
|
|||||||||
|
|
|
× |
1+α(1− μs ) |
|
|
||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||
а энергию после рассеяния |
|
|
|
|
E |
|
|
|
|
|
|
|
|
|
||||
|
|
|
Es = |
|
|
|
|
|
|
|
. |
|
|
|
|
|
||
|
|
|
1+α |
(1− μs ) |
|
|
|
|
|
|||||||||
|
|
|
|
|
|
|
|
|
|
|
||||||||
При |
этом, как и |
ранее, |
ε обозначает энергетический |
эквивалент массы покоя электрона (0.511 МэВ), r0 – классический радиус электрона (2.82×10-13 см), а ne (r ) – концентрацию электронов в точке r . Поскольку в формуле для
151
дифференциального сечения (11.14) присутствует дельта-функция, интеграл (11.13) по E′легко вычисляется:
∞ |
r |
|
|
rr |
|
−rr |
|
r |
e−τ (rr→rrd , Es ) |
|||||||
J&рас = ∫ ∫ ∫ |
ΣKN (Ω → |
|
d |
|
|
| r, E) |
|
|
|
|
|
|
|
× |
||
|
rrd −rr |
|
Σ(rr, E) |
|
rr |
− rr |
|
2 |
||||||||
|
|
|
|
|||||||||||||
∞ 0 4π |
×ψ |
r |
|
r |
|
r |
|
|
d |
|
|
|
|
|||
|
|
|
|
|
|
(11.15) |
||||||||||
|
(r, E, |
Ω) dr dE dΩ |
|
|
|
Обратим внимание, что полученная формула по виду полностью соответствует интегралу (7.1). Под знаком интеграла присутствует
плотность входящих столкновений ψ(rr, E, Ω) , а все остальное
является оценкой, которую надо вычислять после получения точки очередного столкновения:
|
|
r |
r |
r |
|
|
r |
rr |
− rr |
r |
g |
|
(r |
, E,Ω,r |
) = Σ |
|
(Ω → |
rd |
r |
| r, E) |
|
|
R |
|
|
d |
|
KN |
|
rd − r |
Σ(rr, |
E) rrd − rr 2 . (11.16)
Локальная оценка потока (в англоязычной литературе принят термин flux-at-a-point) для вычисления интегральной плотности потока частиц в точке rd , создаваемая за всю историю частицы, представляет собой сумму по всем столкновениям
l |
r |
r |
r |
(11.17) |
εR = ∑gR (rk , Ek , Ωk , rd ) . |
k=1
Отметим, что за возможность вычислять характеристики поля без какого-либо пространственного усреднения приходится платить ухудшением статистических свойств полученной оценки. В самом деле, присутствие в знаменателе формулы (11.16) квадрата расстояния от точки столкновения до точки детектирования в случае нахождения детектора в среде переноса излучения делает оценку неограниченной.
Этот факт, как и в случае оценки по пересечениям (11.10), сам по себе не мешает ее использованию, поскольку математическое ожидание величины оценки конечно. В последнем можно убедиться путем следующих рассуждений. Обозначим через p(rr) плотность вероятности в результате моделирования
столкновений получить случайную точку r . Положим, что эта функция не имеет в окрестности точки детектора rrd каких-то
152
особенностей. Тогда конечность математического ожидания локальной оценки будет определяться интегралом вида
MεR ~ ∫ |
|
1 |
|
|
r r |
, |
|
|
|
|
|
|
p(r)dr |
||
|
|
rr −rr |
|
2 |
|||
|
|
||||||
V |
|
d |
|
|
|
|
взятого по малому объему V вокруг детектора. Этот объем, без ущерба для общности, можно считать сферой малого радиуса δ с центром в точке rrd . Заменяя под знаком интеграла функцию p(r )
ее значением в центре сферы и переходя к сферическим координатам, сразу получим ограниченность математического ожидания:
MεR ~ p(rrd )δ∫4ρπ2 ρ2dρ < ∞.
0
К сожалению, подобный интеграл для дисперсии локальной оценки расходится
δ |
4π |
δ |
||
DεR ~ p(rrd )∫ |
ρ2dρ ~ p(rrd )∫ |
4π |
dρ . |
|
ρ4 |
ρ2 |
|||
0 |
0 |
|
|
При этом в 0 он обладает гиперболической (первой степени) расходимостью, что, несомненно, хуже, чем логарифмическая расходимость оценки по пересечениям.
11.2.3. Другие локальные оценки
Громоздкий вид локальной оценки (11.16) возник в результате того, что рассматриваемая функция детектора содержала сингулярность – пространственную дельта-функцию, поскольку детектор предполагался точечным. В результате, вместо
того, чтобы строить оценку на основе выражения J& =ψ, P* , пришлось применить разложение (5.33) и при n = 1 вычислять компонент рассеянного излучения в виде J&рас =ψ1*,ψ .
Входящий в него компонент ценности столкновений ψ1* никаких
сингулярностей не содержит.
Теперь представим себе детектор, который является не только точечным, но и избирательным по энергии:
153
r |
r |
r |
r |
P(r, E, Ω) = δ (r |
− rd ) δ (E − Ed ) . Оценку для него можно также |
построить с помощью разложения (5.33), но составляющую рассеянного излучения придется считать на основе плотности выходящих столкновений
J&рас = χ2*, χ . |
(11.18) |
В данной формуле присутствует второй компонент ценности выходящих столкновений, который представляет собой следующий интеграл
* r |
r |
∞ |
r |
r′ |
r |
r |
′ r′ |
r′ |
|
||||||||
χ2 (r , E,Ω) = ∫ ∫ ∫∫T (r |
→ r |
| E,Ω)C(E,Ω → E ,Ω |
| r )× |
∞0 4π∞
×T (rr′→ rr′′| E′, Ωr′) P* (rr′′, E′,Ωr′)drr′dE′dΩ′drr′′.
Сучетом конкретного вида функции детектора, а также используя δ-функции, входящие в транспортные ядра, этот 9-
мерный интеграл приводится к одномерному интегралу вдоль луча r
Ω:
* |
r |
r |
∞ exp(−τ0 |
(R, E) −τ |
1(R, Ed )) |
r |
r r |
r |
|||||||
χ2 |
(r , E,Ω) = ∫ |
|
|
|
|
|
|
|
|
|
Σs (E,Ω → Ed ,ϖ R | r |
+ RΩ)dR. |
|||
|
|
r |
r |
r |
|
|
2 |
|
|||||||
|
|
||||||||||||||
|
|
|
0 |
|
|
r |
+ RΩ − rd |
|
|
|
|
|
|
(11.19)
В последней формуле введено обозначение для оптических расстояний
τ0 (R, E) =τ(rr → rr + RΩ, E); τ1(R, Ed ) =τ(rr + RΩ → rrd , Ed ) ,
а также для единичного вектора направления:
ωrR = (rrd − rr − RΩr) / rrd − rr − RΩr (рис.11.3).
Рис.11.3. К построению двойной локальной оценки.
154
На основе интеграла (11.19) получим окончательное выражение для двойной локальной оценки (по пространству и энергии), которое вычисляется каждый раз, когда фазовое состояние частицы соответствует плотности выходящих столкновений, т.е. сразу после розыгрыша новой энергии и направления движения. При этом учтем жесткую связь между углом рассеяния θs и изменением энергии, позволяющую вовсе избавиться от интеграла по лучу:
g |
|
r |
r |
r |
|
) = |
Σ |
s |
(E → E |
d |
| rr+ RΩ) |
× |
|
|
RE |
(r , E,Ω, r , E |
d |
|
|
|
|
|
|||||||
|
|
|
d |
|
|
|
2π ρ sinα sinθs |
|
. |
(11.20) |
||||
|
|
|
|
|
|
|
|
|
|
|||||
|
|
|
|
|
|
|
×exp(−τ0 (R, E) −τ1 (R, Ed )) |
|
|
Входящий в данную формулу угол θs вычисляется из условия, что такое рассеяние переводит энергию E в Ed. Может оказаться, что
при конкретных параметрах rr, E, Ω, rrd , Ed такой угол
отсутствует, тогда и оценку в этой точке не подсчитывают. С учетом сказанного, окончательное выражение для двойной локальной оценки как сумма при разных столкновениях есть
l |
r |
r |
r |
(11.21) |
εRE = ∑gRE (rk , Ek ,Ωk , rd , Ed ) . |
k=1
Дисперсия оценки (11.21) логарифмически расходится по «прицельному параметру» ρ sinα , с которым частица пролетает
мимо точечного детектора.
Действуя аналогично, т.е. используя выражения J&рас =χn* , χ или J&рас =ψn*,ψ нужной (минимальной) степени n,
можно построить оценки и для других типов локальных детекторов. Например, для точечного моноэнергетического мононаправленного детектора. Мы не будем этого делать, поскольку выражение оценки получается громоздким и практического значения не имеет. В главе 13 будет показано, что в подобных случаях более удобным может оказаться моделирование сопряженного уравнения переноса.
В заключение приведем на рис.11.4 общую алгоритмическую схему моделирования фазовых состояний с указанием точек вычисления всех полученных оценок.
155
Оценки по столкновениям εCO и пересечениям εCR , а также
локальная оценка εR вычисляются после розыгрыша транспортного ядра, как только установлена очередная точка столкновения. Оценка по пробегам εPL также вычисляется после
Т-ядра, но может существовать и при вылете частицы за границы системы, т.е. когда фактически точки столкновения нет.
Рис.11.4. Схема алгоритмических точек вычисления различных оценок
Оценка по поглощениям εAB , как уже упоминалось выше,
определяется в процессе розыгрыша C-ядра, если при этом столкновение заканчивается поглощением частицы.
Двойная локальная оценка εRE должна быть рассчитана
перед розыгрышем транспортного ядра. При этом для частиц, только что рожденных в источнике, она даст вклад однократно рассеянного излучения, а в остальных случаях – компонент многократного рассеяния.
Контрольные вопросы
Перечислите основные оценки метода Монте-Карло.
Назовите локальные оценки метода Монте-Карло.
Напишите вид оценки по столкновениям для расчета средней плотности столкновений в объеме детектора.
Назовите достоинства и недостатки локальных оценок.
Каковы статистические свойства локальной оценки потока?
На основе формулы (11.13) напишите локальную оценку
156
потока для столкновения, заканчивающегося образованием электрон-позитронной пары. При этом рассматривайте последующий вылет аннигиляционного излучения как процесс рассеяния первичного фотона (формула (10.7)).
Какие оценки вычисляются многократно в течение истории?
Какие однократно?
157
Глава 12. Методы уменьшения дисперсии
Меж ними все рождало споры И к размышлению влекло…
А.С. Пушкин
§12.1. Неаналоговое моделирование
При рассмотрении в §6.9 вычисления методом Монте-Карло определенных интегралов
J = ∫b f (x) p(x) dx ,
a
подчеркивалось, что такое представление подынтегрального выражения в виде произведения плотности вероятности p(x) и
оценки f(x) не является единственным. Всегда можно перейти к
~
моделированию точек x по другой плотности вероятности p(x) ,
добавив при этом под интеграл дополнительный множитель, равный отношению старой и новой функций:
b |
~ |
p(x) |
|
|
|
||
J = ∫ f (x) p(x) |
|
dx |
|
~ |
|||
a |
|
p(x) |
. |
|
|
При моделировании фазовых состояний частиц в качестве плотности вероятности мы использовали отнормированную должным образом реальную функцию источника (глава 8), или имеющее ясный физический смысл транспортное ядро (глава 9), или же ядро столкновений, в точности отражающее физический процесс столкновения частицы с атомом, ядром или электроном (глава 10). Таким образом, на всех трех этапах плотность вероятности выбираемой фазовой точки соответствовала рассматриваемому физическому процессу. Такое моделирование обычно называют аналоговым.
В отличие от него, идея неаналогового моделирования
предполагает замену реальной (аналоговой) плотности
~
вероятности p(x) на некоторую смещенную плотность p(x) ,
выбранную из математических соображений (уменьшение дисперсии результата, сокращение времени расчета и т.п.). При этом вводится понятие статистического веса частицы, который
158
после каждой выборки нового фазового состояния по смещенной плотности изменяется на величину
= p(x) w ~ r .
p(x)
После n столкновений итоговый вес частицы будет определяться произведением весов, набранных в результате всех произошедших смещенных розыгрышей:
n
Wn = ∏wi .
i=0
Индекс i = 0 присвоен весу, приобретенному в результате смещения функции источника.
Для того, чтобы смоделировать фазовые точки, соответствующие плотности n-х входящих столкновений, следует использовать связь между компонентами неймановского разложения
ψ n (xr) = ∫K (xr′ → xr) ψ n−1 (xr′)dxr′ , |
(12.1) |
которая получена на основе уравнения (5.20). Последовательно применяя ее n-1 раз, будем иметь:
ψ n (xr) = ∫K (xr′ → xr)∫K (xr′′ → xr′)∫K∫K (xr(n−1) → xr(n −2 ) ) ×
(12.2)
× ∫ Q(xr(n ) )T ( xr( n ) → xr(n−1) )dxr′dxr′′Kdxr(n ) .
На последнем шаге плотность первых входящих столкновений расписана (см. формулу (5.21)) через функцию источника и транспортное ядро с использованием их обобщенной записи в фазовых координатах.
Неаналоговое моделирование предполагает следующие смещения.
Для |
розыгрыша |
фазовой точки xr(n) |
рождения частицы |
||
вместо истинной |
функции источника Q(xr(n) ) используются |
||||
замена |
~ r(n) |
) |
и |
компенсирующий |
такое смещение |
Q(x |
статистический вес
= Q(xr(n) ) w0 Q~(xr(n) ) .
159
Для |
получения |
точки xr(n−1) |
первого |
||||||
истинное |
транспортное |
ядро |
T (xr(n) → xr(n−1) ) |
||||||
смещенным |
~ r(n) |
r(n−1) |
) |
и |
вводится |
||||
T (x |
→ x |
|
|||||||
статистическим весом |
|
T (xr(n) → xr(n−1) ) |
|
||||||
|
|
w1 = |
. |
||||||
|
|
~ r(n) |
r(n−1) |
) |
|||||
|
|
|
T (x |
|
→ x |
|
|
столкновения
заменяется
компенсация
Для каждого последующего шага получения точки |
xr(n−i) |
|||||||||||
i-го столкновения |
вместо аналогового |
кинетического |
|
ядра |
||||||||
r(n−i+1) |
r(n−i) |
) |
используется смещенное |
~ r(n−i+1) |
|
r(n−i) |
) |
|||||
K(x |
→ x |
K(x |
→ x |
|
||||||||
с добавлением статистического веса |
|
|
|
|
|
|
|
|||||
|
|
|
wi = |
K (xr(n−i+1) |
→ xr(n−i) ) |
. |
|
|
|
|
|
|
|
|
|
~ r(n−i+1) |
r(n−i) |
) |
|
|
|
|
|
||
|
|
|
|
K (x |
→ x |
|
|
|
|
|
|
Таким образом, к моменту вычисления оценки частица имеет вес, не равный единице. Модифицированные с учетом неаналогового моделирования оценки по столкновениям, пробегам, пересечениям, а также локальная оценка будут иметь вид:
l |
k |
|
ε = ∑εk ∏wi . |
(12.3) |
|
k=1 |
i=0 |
|
В данной формуле εk обозначает конкретное выражение, вычисляемое при каждом столкновении. Оно зависит от вида оценки и искомого функционала (см. 11.1.1, 11.1.3, 11.2.1, 11.2.2).
Буквой l, как и ранее, обозначено последнее столкновение частицы перед обрывом ее истории. Для оценки по поглощениям суммирование отсутствует:
|
|
|
|
|
|
l |
|
|
|
|
|
|
|
|
|
|
ε = εl ∏wi . |
|
|
|
(12.4) |
||
|
|
|
|
|
|
i=0 |
|
|
|
|
|
|
При выборе смещенных ядер или функций плотности |
||||||||||
вероятности |
~ r |
r |
r |
~ r |
r |
r |
r |
~ r |
r |
||
r |
→ |
|
|||||||||
Q(x) |
Q(x) |
; T (x′ → x) →T (x′ → x) ; |
K(x′ → x) → K(x′ → x) |
||||||||
для |
розыгрыша |
фазовой |
точки |
xr |
необходимо |
учитывать |
следующие ограничения. Смещенные ядра и функции плотности вероятности могут обращаться в 0 только в двух случаях:
в этой точке соответствующая аналоговая функция (или ядро)
равна 0;
160