CC BY
88
УДК 539.1
МОДЕЛИРОВАНИЕ ПОТЕРЬ ЭНЕРГИИ ТЯЖЕЛЫХ ЧАСТИЦ С ПОМОЩЬЮ ПАКЕТА ПРОГРАММ Geant4
А. В. Багуля, М. С. Владимиров, В. Н. Иванченко1, Н. И. Старков
В настоящее время Géante является стандартным средством для моделирования экспериментов в физике высоких энергий. В данной работе приводится описание усовершенствования модели процессов ионизации для адро-нов, ионов и классического магнитного монополя, реализованных в Geant4• Обсуждаются результаты моделирования, анализируется их точность в зависимости от выбранных параметров и проводится сравнение с экспериментальными данными.
1. Пакет программ Geant4 [1] позволяет проводить полномасштабное моделирование современных крупных экспериментов в физике высоких энергий (ФВЭ) и астрофизических установок методом Монте-Карло. При разработке Geant4 были использованы современные объектно-ориентированные технологии и язык программирования С++, что является основным отличием от предшествующего пакета программ Geant3. Разработка пакета Geant4 была начата в 1994 году в связи с началом подготовки проекта Большого Адронного Коллайдера (БАК). Первая публичная версия программы была подготовлена в 1998 году. С тех пор Geant4 прошел большой путь совершенствования и параллельно использовался во многих приложениях [2]. Пакет Geant4 являлся основным средством моделирования эксперимента ВаВаг практически с начала набора данных [3]. В настоящее время он принят в том же качестве для будущих эксперименов на БАК, в частности, ATLAS [4] и CMS [5]. Geant4 широко используется в космических [6] и медицинских [7] исследованиях, спектр его применений постоянно увеличивается.
^колого-аналитическая ассоциация "Эко-аналитика", Москва 199992, Россия.
Существенная часть пакета веаг^ - это моделирование физических процессов, организованное как набор программ, предназначенных для моделирования электромагнитных, адронных, оптических и других взаимодействий для очень широкого энергетического спектра [1, 2]. В данной работе описаны модели ионизационных потерь тяжелых частиц, включая моделирование потерь энергии классического магнитного монополя. Показано, что применение сплайн-интерполяции при вычислении потерь энергии, сечений и пробегов позволяет обеспечить большую стабильность и точность результатов. Проведено сравнение предсказаний моделирования с экспериментальными данными.
2. Моделирование потерь энергии в СеаЫ4. Дифференциальное сечение рассеяния протона на свободном электроне среды
Т2
МЕ,Т) 2 а ^ 1
1-/^—+
Тт 2Е2
(1)
¿Т е р2 Т2
где Е, г^ (3 - кинетическая энергия, заряд и скорость налетающей частицы; Т - кинетическая энергия дельта-электрона; ге, тпе - классический радиус и масса электрона; Тт - максимальная энергия, которая может быть передана электрону,
Тт = 2те(72 - 1)/(1 + 27(те/М) + (те/М)2), (2)
где 7 - релятивистский фактор налетающей частицы, М - ее масса. Сечение образования дельта-электрона можно записать как
<т(Е,Тс)=Т/^<1Т, (3)
тс
где Тс - значение порога образования дельт а-электронов. Введение порога необходимо как для оптимизации скорости моделирования, так и для обеспечения условий, при кото рых выражение (1) справедливо (приближение квазисвободного электрона). Тормозная способность (удельные ионизационные потери энергии) для данной частицы может быть представлена в виде
Ттп ,
&(£?) = ад - / Т^с/Т, (4)
Тс
где - полные потери энергии, а 3Г(Е) - ограниченные потери энергии (переда-
ча энергии меньше Тс). При реальном моделировании используются выражения (3) и (4), причем с уменьшением величины порога Тс возрастает сечение рождения дельта-электронов и уменьшается величина средних потерь энергии на шаге частицы.
Рассмотрим как реализовано моделирование потерь энергии ионов в пакете Geant4. Для ускорения расчетов перед началом моделирования создаются таблицы потерь энергии 5(Т), пробегов R{T), обратных пробегов T(R) и сечений сг(Т). В целях экономии оперативной памяти компьютеров и времени для построения таблиц в случае положительно заряженных адронов таблицы потерь энергии и пробегов создаются только для протонов и альфа-частиц. Эти таблицы используются в дальнейшем для вычисления потерь энергии более тяжелых заряженных частиц. Такой подход принципиально важен, поскольку количество стабильных изотопов превосходит 3000 и при взаимодействии адронов с ядрами может образоваться любой из них.
Моделирование в пакете Geant4 осуществляется путем розыгрыша методом Монте-Карло дискретных шагов частицы в веществе. Вычисление средних потерь энергии АЕ после любого шага производится следующим образом:
AE = Srx, (Srx<Et)
AE = E-E{R{E)-x), (Srx>E(),
где x - длина шага, £ - параметр, определяющий условие использования линейного приближения.
Вследствие вероятностного характера взаимодействия частицы со средой потери энергии на данном отрезке траектории испытывают заметные флуктуации. В Geant4 для описания флуктуаций применяются две различные модели, зависящие от параметра который равняется нижнему пределу среднего числа взаимодействий частицы на шаге моделирования. В случае достаточно большого числа взаимодействий
АЕ > кТс и Тт< 2Тс (6)
флуктуации распределены по нормальному закону с дисперсией
П2 = ((АЕ - (АЕ))2} = 2Trlmce2N^Tcx il - Çj , (7)
где Ne\ - электронная плотность среды. Если условия (6) не выполняются, то применяется более сложная модель флуктуации потерь энергий в тонких поглотителях [8]. В зависимости от числа взаимодействий на шаге реализуются распределения Вавилова или Ландау.
2.1. Потери энергии адронов. Адроны, проходя через вещество, теряют энергию в основном на ионизацию. Среднее значение ограниченных потерь энергии (4) в случае
достаточно быстрых тяжелых частиц задается при помощи формулы Бете-Блоха [9] с поправками
+ /Ъ, (8)
где I - средний потенциал ионизации атома, 6 - поправка, учитывающая эффект плотности, Се/£ - член, учитывающий взаимодействие с электронами на внутренних оболочках атомов, Z- атомный номер, Ен - поправки высших порядков.
Для скорости заряженного адрона (3 < 0.05 (энергия 1 МэВ для протона) формула (8) становится неточной. В этом случае скорость адрона становится сравнимой со скоростью атомных электронов. Поэтому для протонов формула (8) применяется при Т > 2 МэВ. При меньших энергиях в Сеап14 применяется параметризация потерь энергии. Для 72 материалов, доступных в базе данных Национального института по стандартам и технологиям США (ШБТ) [10], прямо используются ШБТ таблицы. Для остальных материалов применяется параметризация [11], предложенная Международной комиссией по радиологическим единицам и измерениям (1С1Ш):
= Е < 10 кэВ,
с _ ^оуу^Ы^ с _ л т-|0.45
^е — т; 5 — ,
"г ^Ъhigh
5ы6ь = (1 + ^ + , 10 кэВ < Е < 2 МэВ, (9)
где б^е _ полная тормозная способность, связанная с передачей энергии электронам, Аг - пять параметров, используемых для аппроксимации, причем они индивидуальны для протонов и а-частиц для каждого атома в диапазоне атомных чисел от 1 до 92. Отметим, что при очень низких энергиях (/? < 0.01) модель свободного электронного газа [12] предсказывает, что торхмозная способность пропорциональна скорости адрона, что согласуется с используемой параметризацией (9).
Для низких энергий учитываются также неионизационные (ядерные) потери 5П, параметризация которых проводится согласно рекомендации 1С1Ш [13]. Полные энергетические потери
St = Se + Sn. (10)
Точность формулы Бете-Блоха оценивается в 2-3% для энергий Е > 5 МэВ, точность вычислений потерь энергии протонов в интервале энергий 100 кэВ < Е < 5 МэВ оценивается в 5%, при более низких энергиях ухудшается с 5 до 20% при 1 кэВ.
2.2. Потери энергии ионов. Для не очень больших энергий 7те/М << 1 сечение (3), ионизационные потери (8) и (9) не зависят от массы адрона, а только от скорости. Исключение составляет поправка к формуле Бете-Блоха (8). Поэтому потери энергии 5в1 (Т) заряженного адрона или иона с кинетической энергией Т могут быть выражены через потери энергии протона, движущегося с той же самой скоростью
5е1(Т) = (Тр) + ЩТ) - Гр(Тр)) Тр = Т(МР/М), (11)
где Zetf - эффективный заряд иона, - масса протона, 5ер - тормозная способность протона, которая соответствует "масштабированной" кинетической энергии Гр. Величины Зер(Т) и пробеги вычисляются из таблиц, поправки высокого порядка /', вычисляются на каждом шаге моделирования.
Также на каждом шаге вычисляется эффективный заряд иона. При высоких скоростях ион целиком лишен электронов, поэтому его заряд не отличается от заряда ядра. Для медленных ионов ситуация кардинально изменяется, потому что начинается захват электронов атомов вещества на электронные оболочки иона, что приводит к уменьшению его заряда, в результате чего его потери энергии уменьшаются. При движении в плотной среде наряду с захватом электронов происходит и передача электронов атомам среды, что приводит, с одной стороны, к установке равновесия между ионом и средой, с другой - к дополнительным флуктуациям потерь энергии из-за флуктуаций заряда иона [14]. До настоящего времени отсутствует достаточно развитая теория этого процесса, но предложено несколько количественных эмпирических соотношений. Эффективный заряд иона выражается через обычный заряд иона Zeя = 7а величина 7г параметризована в соответствии с формулами, предложенными Зиглером [13].
2.3. Потери энергии магнитного монополя. С момента предсказания магнитного монополя предпринимались многочисленные попытки экспериментально его зарегистрировать. Тяжёлый магнитный монополь должен обладать высокой проникающей способностью и создавать на своём пути сильную ионизацию. Это один из признаков, по которому его можно зарегистрировать. Большинство теоретических предположений предсказывают массу магнитного монополя на уровне ГэВ. Но некоторые тео рии свидетельствуют в пользу существования магнитного монополя с гораздо меньшей
массой. В лаборатории им. Ферми (США), где проводится эксперимент по поиску магнитных монополей, удалось установить новый нижний предел на их массу в 360 ГэВ [15]. Предполагается, что поиски монополя продолжатся в экспериментах на БАК. Поэтому естественным было желание добавить в пакет веап14 описание процессов ионизации, производимой монополем.
Тс = 350 кэВ Гс-1.7каВ
80 90 100 Ц мм
Рис. 1. Сравнение результатов моделирования с помощью пакета СеапЦ (версия 9.1,) распределений энерговыделения при прохождении пучка протонов с начальной энергией 70 МэВ вдоль водяной мишени для различных значений порога образования дельта-электронов Тс.
Рис. 2. Сравнение результатов моделирования с помощью пакета СеапЦ (версия 9А) распределений энерговыделения при прохождении пучка протонов с начальной энергией 110 МэВ вдоль водяной мишени для различных значений
Магнитный монополь при прохождении через вещество теряет энергию в результате взаимодействия его заряда с электрическим полем атомных электронов. Ионизационные потери монополя вычисляются по формуле С. Алена [16]
¿Е 4тгЛГе1еУ
¿х тес
(12)
Эта формула верна при (3 > 0.1. При /3 < 0.01 для всех материалов используется следующая аппроксимация [17]:
-¿Е/±г = (45(ГэВ • см2)/г)02/?. (13)
Для промежуточного интервала скоростей 0.01 < ¡3 < 0.1 в настоящий момент реализована линейная интерполяция.
Рис. 3. Вычисление пробегов по таблицам для случая линейной интерполяции (случай малого шага моделирования).
Рис. 4. Распределение энерговыделения при прохождении пучка электронов с начальной энергией 1 МэВ вдоль алюминиевой мишени, промоделированное с помощью пакета СеаЫУ, (версия 9.1) при отключении процессов многократного рассеяния и флуктуации.
3. Оптимизация процесса моделирования радиоактивной дозы от пучка тяжелых частиц. Точность результатов веаМ4 и скорость расчетов зависят от значений параметров моделирования, поэтому актуальным становится их правильный выбор. В данной работе изучалось влияние параметров моделирования на форму распределения потерь энергии протонов и электронов. Использовалось приложение Нас1г01 пакета веап14 [18], разработанное авторами данной работы в качестве примера расчетов замедления пуч ка тяжелых частиц в протяженной мишени. Мишенью служил цилиндр, разделенный на тонкие слои перпендикулярно к направлению пучка. Потери энергии усреднялись в каждом слое. Размеры цилиндра и число слоев могут быть заданы пользователем. В качестве материала мишени были выбраны вода, алюминий и свинец.
Когда в качестве первичной частицы выбирался протон, моделирование проводи лось при энергиях 70, 110, 160 и 400 МэВ. Сначала определялось влияние величины Тс на форму кривой энерговыделения (рис. 1). При меньших значениях Тс происходило сглаживание Брэгговского пика. Это означает, что в случае явного моделирования
дельта-электронов наблюдается больше флуктуации. Возможно, это объясняется тем, что в этом случае моделирование флуктуаций потерь энергии не является корректным. Было установлено, что результаты различаются для малых и больших значений параметра Тс, причем для меньших значений энергии это различие проявляется больше.
Рис. 5. Сравнение результатов моделирования с помощью пакета СеапЬ^ (версия 9.26е/а) энерговыделения пучка протонов с начальной энергией 100 МэВ (а) и ионов углерода с начальной энергией 100 МэВ/а.е.м. (б) в водяной мишени с экспериментальными данными [20] (сплайн-опция активирована).
Когда шаг моделирования достаточно мал (0.1 мм в нашем случае), при малых значениях £ (см. (3)) на распределении потерь энергии по глубине появляются особенности - кривая принимает волнообразный характер (рис. 2). Также в зависимости от £ меняется положение Брэгговского пика. Моделирование показало, что оптимальные значения £ лежат в пределах 0.01-0.001. Этот эффект возникает, когда длина шага х меньше расстояния между узлами таблицы потерь энергии и пробегов (рис. 3). В этом случае применение линейной интерполяции при вычислении потерь энергии (пробегов, сечений) между узлами таблицы приводит к ошибкам. Соответственно, существенное увеличение числа ячеек в таблице (до 1200 против стандартных 120) устраняет численную нестабильность. При этом положение пика становится нечувствительным к параметру £. Но на практике существенное увеличение числа ячеек таблицы не является оптимальным. Увеличение размеров таблиц, которые строятся для всех материалов и основных частиц, приводит к значительному росту занимаемой приложением оперативной памяти
Рис. 6. Сравнение рассчитанных с помощью пакета СеапЦ (версия 9.1,) зависимостей ионизационных потерь от энергии для магнитных монополей массой 100 и 500 ГэВ (а) и протона (б).
Е, МэВ
Рис. 7. Сравнение рассчитанных с помощью пакета СеапХ^ (версия 9А) зависимостей пробега от энергии для протона и магнитного монополя.
компьютера. Поэтому такое решение возможно только для относительно простых приложений. Поэтому в веап14 был реализован метод интерполяции физических таблиц на основе кубического сплайна [19]. В результате применения этого метода удалось достичь стабильности моделирования положения Брэгговского пика, а также стабильности его формы при существенной вариации порогов образования дельта-электронов.
600 ^ 400 200 0
Еп = 50 ГэВ
— А1
......жидкий Аг
11. ■ 11
1 ■ ■ * ■1 ■ ■ ■
\
I I I М I I I I I I I I I I I I I I Л I I I I I
О 10 20 30 40 50 60 70
80 90 100 мм
120 140
I., мм
£0 = 300 ГэВ -А1
.......жидкий Аг
150 200 250
300 350 400 мм
2400
2200 2000 21800 _11
£0= 1000 ГэВ — А1
21600
$1400
^ 1200
¿ЕМ о о с
400
200 !
0 ........ ■ 1 , 1 . . . 1 1 ^^
200
400 600
800 1000 Ь, мм
Рис. 8. Сравнение результатов моделирования с помощью пакета СеапЦ (версия 9.1) распределений энерговыделения магнитного монополя массой 100 ГэВ в алюминии и жидком аргоне при различной начальной энергии Е0) полученных после применения интерполяции.
Начиная с версии 9.2Ье1а, пользователь Сеап14 может выбирать, какой метод интерполяции использовать. Причем стоит отметить, что сплайн-интерполяция физических таблиц реализована в Сеап14 таким образом, что ее использование не приводит ни к увеличению времени моделирования, ни к дополнительному расходованию оперативной памяти компьютера.
Такие же исследования проводились и для электронов при энергии 1 МэВ. Выбор энергии определялся тем, что, во-первых, электроны с этими энергиями составляют значительную часть электромагнитных ливней, а во-вторых, на эти энергии приходится максимум рентгеновского тормозного излучения. Аналогичный эффект наблюдался при отключении процессов флуктуации и многократного рассеяния, которое все сглаживает
(рис. 4). Но для точного моделирования потерь энергии в тонких поглотителях этот эффект следует учитывать. Применение сплайн-интерполяции для электронов также представляется полезным.
4. Результаты. Сравнение результатов моделирования Брэгговского пика в сравнении с экспериментальными данными для различных ионных пучков очень важно для ФВЭ, так как заметная часть адронных ливней выделяется вследствие процессов ионизации адронов/ионов. Это существенно и для космических, и особенно для медицинских приложений, так как оптимизация радиационной терапии является одним из основных современных направлений исследований в физической медицине. На рис. 5 представлены результаты моделирования прохождения углеродного пучка с энергией 100 МэВ/а.е.м. и протонного пучка с энергией 110 МэВ через водяную мишень. Результаты моделирования хорошо воспроизводят экспериментальное распределение дозы по глубине. Использование модели парного (Binary) каскада [20] позволяет хорошо описать фрагментационный хвост распределения энерговыделения за Брэгговским пиком.
С целью изучения особенностей торможения монополей в различных средах было разработано специальное приложение "monopole", которое включено в пакет Geant4 в качестве примера. С помощью данного приложения получены графики потерь энергии и пробегов монополей в алюминии и жидком аргоне (рис. 6-8), являющимися основными материалами электромагнитного калориметра детектора ATLAS.
5. Заключение. В данной работе описаны методы моделирования потерь энергии протонов, ионов и магнитных монополей, реализованные в рамках пакета программ Geant4. Разработаны методы оптимизации использования таблиц энергетических по терь и пробегов. Учтены поправки к основным формулам энергетических потерь. В результате достигнута устойчивость результатов моделирования положения Брэгговского пика, результаты хорошо согласуются с экспериментальными данными. Подготовленное моделирование потерь энергий классическим магнитным монополем позволит провести поиск тяжелых монополей в экспериментах на БАК.
Статья подготовлена при частичной поддержке РФФИ, грант 09-02-91065.
ЛИТЕРАТУРА [1] S. Agostinelli et al. (GEANT4 Collaboration), Nucí. Instr. and Meth. A 506, 250 (2003).
[2] J. Allison et al., IEEE Trans. Nucl. Sei. 53, 270 (2006).
[3] S. Benerjee et al., J. Phys. Conf. Ser. 119, 032007 (2008).
[4] B. Dowler et al., Nucl. Instr. Meth. A 482, 94 (2002).
[5] P. Arce et al., Nucl. Instr. Meth. A 502, 687 (2003).
[6] G. Santin, V. Ivanchenko, H. Evans et al., IEEE Trans. Nucl. Sei. 52, 2294 (2005).
[7] S. Jan et al., Phys. Med. Biol. 49, 4543 (2004).
[8] K. Lassila-Perini, L. Urban, Nucl. Instr. Meth. A 362, 416 (1995).
[9] C. Caso et al., Europ. Phys. Jour. С 3, 1 (1998).
[10] http://physics.nist.gov/PhysRefData/
[11] A. Allisy et al., ICRU Report 49, 1 (1993).
[12] J. Linhard and A. Winther, Mat.-Fys. Medd. 34, 4 (1964).
[13] J. F. Ziegler and J. M. Manoyan. Nucl. Instr. and Meth. В 35, 215 (1988).
[14] M. А. Кумахов, Ф. Ф. Комаров, Энергетические потери ионов в твердых телах (Минск, Изд-во БГУ, 1979).
[15] A. Abulencia, Phys. Rev. Lett. 96, 201801 (2006).
[16] S. P. Ahlen, Rev. Mod. Phys. 52, 121 (1980).
[17] G. R. Kalbfleisch et al., Phys. Rev. D 69, 052002 (2004).
[18] A. Bagulya, I. Gudowska, V. Ivanchenko, N. Starkov, Proton/ion Bragg peak validation. 11th Geant\ Collaboration Workshop (Lisbon, 9-14 October 2006). http: / / geant4. web.cern.ch/geant4/collaboration / working, groups/electromagnetic/ presentations.shtml
[19] W. H. Press et al., Numerical recipes in C: the art of scientific computing (Cambridge University Press, 1997).
[20] G. Folger, V. N. Ivanchenko and H. P. Wellisch, Eur. Phys. J. A 21, 407 (2004).
[21] U. Oelfke, G. K. Y. Lam, M. S. Atkins, Phys. Med. Biol. 41, 177 (1996).
Поступила в редакцию 20 октября 2008 г.
