Занятие 3. Задание физических процессов в Geant4.

На мой взгляд, это самое сложное, с чем может столкнуться пользователь при моделировании своего эксперимента. Тем не менее, я постараюсь изложить данный этап моделирования как можно более понятно и по возможности подробно. Повторюсь, что если у читателя возникнут вопросы, я всегда готова ответить на них в комментариях.

Как и в предыдущем посте про физические модели, сегодня мы будем говорить о физике применительно к адронной терапии. Примерный набор физических процессов,  которые используются в этой области физики можно посмотреть в примере "hadrontherapy" в директории "/examples/advanced/hadrontherapy" внутри дистрибутива Geant4 10.01 (файл HadrontherapyPhysicsList.сс). Или в нашем примере - там, по сути, используется та же физика.

Теперь поговорим подробно о способах, которыми эта физика задается. Всю систему моделирования физики внутри Geant4 можно разделить на три принципиальных элемента:

1. Физическая модель - генерирует конечное состояние частицы после взаимодействия и обновляет параметры моделируемой частицы.
2. Физический процесс - содержит необходимые таблицы поперечных сечений и используемые физические модели.
3. Список физических процессов (Physics List) - как понятно из названия, элемент, который содержит в себе набор физических процессов для каждой частицы.

При написании программы пользователь обязательно должен создать свой набор физических процессов в своем классе, который должен быть унаследован от G4VUserPhysicsList. После чего указатель на него он передает в G4RunManager. Вот как это должно выглядеть в коде в функции main:

G4VUserPhysicsList* physicsList = new PolyethylenePhysicsList;
runManager->SetUserInitialization(physicsList);

Но есть несколько серьезных упрощений. Во первых, Geant4 предоставляет уже готовые наборы физики (Reference Physics List), предназначенные для различных областей. Например нам более менее подходит QGSP_BERT, ориентированный на эксперименты с адронами. Какие модели туда входят и, вообще какие готовые наборы физики нам подходят вы сможете понять отсюда. Обратите внимание на раздел "Medical and industrial neutron applications" и "Low energy dosimetric applications". Вообще нам подходят разделы с ключевыми словами hadron и low energy.
Такие готовые наборы физики подключаются очень просто. Вот так:

runManager->SetUserInitialization(new QGSP_BIC);

Только не забудьте сделать include заголовка класса, в данном случае QGSP_BIC.hh.

Второй по убыванию сложности метод задания физики - это наследование класса G4VModularPhysicslist, который, в свою очередь является наследником класса G4VUserPhysicsList. У нас получается две ступени наследования!
В функции main это будет выглядеть так:

  G4VModularPhysicsList* physicsList = new PolyethylenePhysicsList;

Класс G4VModularPhysicsList создан для того, чтобы мы могли построить физику из набора блоков, о которых мы говорили, когда обсуждали физику в предыдущем посте: G4EmStandardPhysics, G4EmExtraPhysic и так далее. Если помните, то мы выяснили, что каждый из этих классов отвечает за свою область физики. Более того, каждый из этих классов является наследником одного и того же материнского класса G4VPhysicsConstructor. Этот класс содержит две важные функции, которые должен определить по своему каждый класс, который его наследует: ConstructParticle() и ConstructProcess(). Это означает, что обе эти функции по-своему определены в G4EmStandardPhysics и во всех других таких классах. В первой функции создаются все частицы, участвующие в этом разделе физики (в том числе теоретически возможные вторичные), во второй функции к каждой такой созданной частице подсоединяются соответствующие ей физические процессы. Пожалуйста, запомните это и не спутайте дальше с точно так же называющимися функциями в классе G4VUserPhysicsList! О которых мы еще поговорим.

Итак, что у нас есть: класс G4VModularPhysicsList и набор дочерних классов от G4VPhysicsConstructor. Для использования их в нашем G4VModularPhysicsList, в классе G4VModularPhysicsList есть специальная функция RegisterPhysics, в которую как аргумент передается указатель на объект нужного нам класса, наследующего G4VPhysicsConstructor. Взгляните в пост от занятия 2 и 3/4. Там как раз показан список из таких функций RegisterPhysics(). (Обратите внимание, что этот список должен содержаться в конструкторе нашего класса - наследника G4VModularPhysicsList.) Внутри функции RegisterPhysics() переданный в нее класс физики добавляется в массив и хранится там. Что же происходит дальше? Когда мы запускаем расчет, G4RunManager обращается к нашему дочернему классу G4VModularPhysicsList и вызывает в нем две функции по очереди ConstructParticle() и ConstructPhysics(), которые (по доброте душевной разработчиков) называются точно так же, как и функции класса G4VPhysicsConstructor. В чем отличие? В том что эти две функции, пустые, с точно такими же названиями содержатся в классе G4VUserPhysicsList. Когда класс G4VModularPhysicsList наследует класс G4VUserPhysicsList, эти функции с точно такими же названиями переходят в G4VModularPhysicsList. Но здесь они уже не пусты! Класс G4VModularPhysicsList их определяет. В своей функции ConstructParticle() он обращается к массиву физических классов - наследников G4PhysicsConstructor и вызывает у каждого по очереди его собственную функцию ConstructParticle(). То же самое происходит и в функции ConstructPhysics() класса G4VModularPhysicsList - в ней идет обращение к тому же массиву, и на каждый класс, содержащийся там, вызывается его функция ConstructPhysics().

Теперь займемся классом G4VUserPhysicsList и самым сложным способом задания физики. Он сложен потому, что тут мы не используем модули, ту мы вручную будем создавать каждую частицу и к ней присоединять все нужные физические процессы. Вообще, мы можем пойти еще дальше, и к каждому процессу вручную подсоединить нужные физические модели. Особенно это касается адронной физики, где моделей, описывающей ядерные взаимодействия великое множество для каждого диапазона энергий.
Как я уже говорила, в классе G4VUserPhysicsList есть две пустые функции ConstructProcess() и ConstructParticle(), которые нам обязательно надо задать. В предыдущем случае за нас это делал G4VModularPhysicsList, но теперь мы сыграем в Чака Норриса и сделаем все сами.
Сначала один не очень очевидный факт. Задавать упомянутые функции мы будем точно так же, как они задаются в классах типа G4EmStandardPhysics (помните, что у их материнского класса точно такие же по названию функции?) с той лишь разницей, что теперь мы создаем все возможные в нашем эксперименте частицы и подключать к ним все возможные для них физические процессы.
Давайте заглянем в файл-исходник класса G4EmStandardPhysics вот здесь и посмотрим как это делается. С функцией ConstructParticle() все довольно просто. Создаются классы - синглетоны для каждого типа частицы, описывающие ее свойства. Во всей программе возможен только один объект каждого такого класса.
Что же происходит в функции ConstructProcess()? Обратите внимание на следующий кусок кода:

theParticleIterator->reset();
while( (*theParticleIterator)() ){
     G4ParticleDefinition* particle = theParticleIterator->value();
     G4String particleName = particle->GetParticleName();
     if(verbose > 1)
       G4cout << "### " << GetPhysicsName() << " instantiates for " 
              << particleName << G4endl;
 
     if (particleName == "gamma") {
       ph->RegisterProcess(new G4PhotoElectricEffect(), particle);
       ph->RegisterProcess(new G4ComptonScattering(), particle);
       ph->RegisterProcess(new G4GammaConversion(), particle);...

Здесь происходит следующее. Используется объект theParticleInterator, который является static и доступен из любой области программы. Он пробегает по всем инициализированным нами частицам и по очереди возвращает указатель на каждую из них. Это значит, что цикл while выполняется, пока мы не переберем все частицы.

Следующей строчкой после инициализации цикла мы получаем указатель на частицу - particle и в зависимости от того как она называется (функция GetParticleName()), мы подключаем к ней процессы, относящиеся к данному типу частицы. Делаем это мы с помощью множественных if циклов. Обратите внимание сколько их много. Первый из них, для гамма кванта, показан здесь.

Для того, чтобы понять, что происходит в цикле if также потребуются некоторые усилия. Обратите внимание на функцию RegisterProcess! Здесь мы имеем дело с классом G4PhysicsListHelper, а именно его объектом ph. Рядом с ним не стоит название класса, потому что в начале функции ConstructProcess, в которой мы сейчас находимся, указатель на него уже получен (найдите это место). Обратите внимание, что класс G4PhysicsHelper static.

Этот класс существует исключительно для того, чтобы облегчить процесс поключения процессов к частице. Внутри функции RegisterProcess он имеет дело с классом G4ProcessManager самой частицы. Он проверяет какие процессы уже к ней поключены, и поключает новый процесс, ставя ему в соответствие очередь, в которую он вызывается.

Эта очередь - нетривиальная вещь, так как она продиктована физикой процесса. Есть три главных этапа вызова процессов во время одного шага частицы - AlongStepDoIt, PostStepDoIt, AtRestDoIt.  В зависимости от такого какой процесс - непрерывный или дискретный - он добавляется в список к одному (иногда к нескольким) из этих этапов. Например, поглощение частицы в конце трека - это чистый AtRestDoIt процесс, а неупругое столкновение протона с ядром с рождением вторичных частиц - PostStepDoIt процесс. Внутри каждой этой группы есть свой список последовательности вызова процессов.  Все это можно подклучать вручную, но тогда нам не надо использовать G4PhysiсsListHelper, а получать от каждой частицы ее G4ProcessManager. Если вам интересно как это делается и как правильно выстраивать процессы в очередь - об этом написана целая глава в User's Guide for Application Developers  для Geant4 (она доступна в pdf бесплатно здесь). Страница 177, раздел 5.2.

Вернемся к нашему if циклу.  Внутри него мы поключаем к каждому типу частицы процессы ей соответствующие с помощью функции RegisterProcess, в которую в качестве аргументов передаются указатели на частицу и на процесс для нее.

But we need to go even deeper. Каждый процесс имеет внутри себя модель или набор моделей, которые, собственно, всю работу и выполняют. Несколько моделей бывает нужно, если каждая из них работает, например, для своего диапазона энергии. Или если одна модель теоретическая, другая использует экспериментальные таблицы сечений.

Помните, на лекции я не смогла ответить на вопрос чем функция электромагнитного процесса SetEmModel отличасется от AddEmModel? После прочтения километра кода отвечаю. Если вы хотите использовать одну модель, то поключайте ее через функцию SetEmModel обязательно с номером 1 как второй аргумент функции! Если будет 0 или 2 или больше, программа эту модель проигнорирует! AddEmModel используйте если вам надо подключить несколько моделей. Причем цифра второго аргумента должна соответствовать узости специализации каждой модели. Например, если у вас есть две модели, одна из которых для низких энергий, а вторая для широкого диапазона - то модель низких энергий должна идти под цифрой 2, а общая модель под цифрой 1. И не используйте нулей, код программы игнорирует модели под нулевым номером. Странно, согласна.

В завершение я хочу упомянуть, что у электромагнитных процессов есть механизмы ограничения шага частицы, которые повышают точность моделирования. Они предусмотрены и для торможения частицы, и для многократного рассеяния, и общие. Вообще, физику можно использовать и без настройки этих параметров, так как они разумно выставлены по умолчанию. Однако, в специфических экспериментах это может понадобится. По поводу этих алгоритмов я напишу отдельный пост, но, к сожалению, на английском языке. Дело в том, что это может быть интересно иностранному сообществу пользователей тоже.

Занятие 2 и 3/4. Behind the scenes: Обзор физических моделей для протонов, содержащихся в Geant4

Сегодня я хочу изложить подготовительную информацию для нашего следующего, третьего занятия. На третьем занятии мы будем говорить о том, как задавать набор физических процессов и соответствующих моделей в Geant4. Сегодня же я хочу отвлечься от кодинга и немного обсудить физику, заложенную в эти модели, а так же логику их использования при описании эксперимента.
Сначала загляните в файл PolyethylenePhysicsList.cc и обратите внимание на эти строки:

RegisterPhysics(new G4EmStandardPhysics_option3);
RegisterPhysics(new G4HadronPhysicsINCLXX);
RegisterPhysics(new G4EmExtraPhysics);
RegisterPhysics(new G4HadronElasticPhysics);
RegisterPhysics(new G4StoppingPhysics);
RegisterPhysics(new G4IonBinaryCascadePhysics);
RegisterPhysics(new G4NeutronTrackingCut);
//    RegisterPhysics(new G4DecayPhysics);
//    RegisterPhysics(new G4RadioactiveDecayPhysics);

Вы видите, что все применяемые физические модели разбиты на крупные разделы. Среди них есть основные:

• Электромагнитная физика (G4EmStandardPhysics_option3). Она включает в себя основные электромагнитные процессы для протонов: многократное рассеяние на электронных оболочках атомов и ионизационное торможение в веществе. Кроме того, в этот раздел входят все процессы и для других частиц, которые могут родиться в эксперименте как вторичные. Например, процессы для фотонов включают в себя: образование пар (здесь оно немного непривычно для меня называется G4GammaConversion), эффект Комптона, Релеевское рассеяние, фотоэлектрический эффект. И так далее. Каждому типу частицы в этом разделе ставятся в соответствие основные электромагнитные процессы ей соответствующие. Как это все подключается мы обсудим на следующем, третьем, занятии, посвященном программированию.
• Упругое рассеяние адронов (G4HadronElasticPhysics). В этом разделе описывается процесс упругого рассеяния адронов на ядрах вещества, а не на электронных оболочках.
• Неупругое рассеяние адронов (G4HadronPhysicsINCLXX). В этом разделе описывается процесс неупругого рассеяния адронов на ядрах вещества.

Оставшиеся разделы в списке не настолько важны для нас, но, тем не менее, рассмотрим и их:

• Дополнительная электромагнитная физика (G4EmExtraPhysics). В нее входит три подраздела:

1. Синхротронное излучение (G4SynchrotronRadiation). Электромагнитное излучение, испускаемое ультрарелятивистской заряженной частицей, ускоряемой радиально (a ⊥ v). По названию понятно, что такое излучение было впервые обнаружено в ускорителе сихротроне. В нашем случае сихротронное излучение возможно, если быстрые электроны будут двигаться через магнитное поле (мы пробовали его включать на прошлом занятии). Для более тяжелых частиц синхротронное излучение в нашем эксперименте невозможно, так как они не достигают релятивистских энергий (напомню, что энергия протонов в пучке в нашем эксперименте - 160 МэВ).
2. Фотоядерный эффект (G4PhotoNuclearProcess) или поглощение гамма-кванта ядром. При этом ядро получает избыток энергии без изменения своего нуклонного состава. Если переданная ядру энергия превосходит энергию связи нуклона в ядре, то распад образовавшегося составного ядра происходит чаще всего с испусканием нуклонов, в основном нейтронов.
3. Взаимодействие мюона с ядром (G4MuonNuclearProcess) путем обмена виртуальным фотоном.

• Ядерный захват остановившихся отрицательно заряженных частиц (G4StoppingPhysics). Применим для остановившихся  анти-протона, μ-, K-, Σ-, Ξ-, Ω- , анти - Σ+.
• Ядерные рекакции распространенных ионов с ядром (G4IonBinaryCascadePhysics). Моделируется с помощью бинарного внутриядерного каскада. Ионы могут возникать как продукты неупругих взаимодействий протонов с ядрами вещества.
• Ограничитель на моделирование трека нейтронов (G4NeutronTrackingCut).  Так как нейтрон - нейтральная частица, он не взаимодействует с электронными оболочками атомов. Замедление нейтронов происходит при упругих столкновениях с ядрами. Если до столкновения ядро покоилось, то после столкновения оно приходит в движение, получая от нейтрона некоторую энергию, вследствие чего нейтрон замедляется. Однако это замедление нейтронов не может привести к их полной остановке из-за теплового движения ядер. В итоге, необходим ограничитель на трекинг нейтронов по времени или по энергии нейтрона.

Также можно подключить оставшиеся два раздела физики связанные с распадом нестабильных частиц (G4DecayPhysics) и радиактивным распадом ядер (G4RadioactiveDecayPhysics). Однако, эти процессы вносят чрезвычайно малый вклад в результат нашего эксперимента и их можно не учитывать.

Вернемся к основным процессам взаимодействия протонов с веществом и рассмотрим их более подробно.

• G4EmStandardPhysics_option3

При подключении раздела электромагнитной физики, протонам (и антипротонам) приписываются набор определенных физических процессов взаимодействия. Каждый такой физический процесс содержит в себе модель, которая описывает данное взаимодействие. Важно, что один и тот же физический процесс могут описывать разные модели, и эти модели взаимозаменяемы. То есть мы можем поменять модель в этом процессе или подключить несколько моделей, определив для каждой диапазон энергий частиц, в котором она применима.

Итак, обсудим процессы, которые у нас используются в электромагнитной физике:

Упругое рассеяние протонов на электронных оболочках атомов.  Как раз в этом месте есть небольшая сложность. Для этого взаимодействия предусмотрена два варианта процесса, а не модели. Дело в том, что Geant4 позволяет описывать это взаимодействие двумя разными способами.

Первый - это детализированное описание рассеяния в отдельности. Такой процесс так и называется G4CoulombScattering. К нему подключается модель G4eCoulombScatteringModel. В этом процессе, каждое упругое рассеяние моделируется отдельно. Преимущество этого процесса в том, что он использует таблицы сечений взаимодействия, полученные теоретически, в отличие от второго варианта, процесса многократного рассеяния. Недостатком является то, что из-за того что моделируется каждое взаимодействие, количество шагов (G4Step) частицы сильно возрастает по сравнению с многократным рассеянием. Следовательно, возрастает и время расчета, в сотни, тысячи раз. Такой процесс имеет смысл подключать только если мы имеем дело со средой низкой плотности, например, воздухом. Физическая модель G4hCoulombScatteringModel основана на теории рассеяния Вентзеля (Wentzel). Подробное описание теории можно посмотреть в Geant4 Physics Reference Manual стр. 67.

Второй способ описания - многократное рассеяние, процесс G4hMultipleScattering (напомню, что буква h в названиях процессов, означает что они применимы к адронам). К нему могут быть подключены две модели, каждую из которых мы рассмотрим отдельно:

1. G4UrbanMscModel - модель Урбана (Urban), которая применима для всех типов частиц. Модель Урбана использует специально подобранные функции углового и пространственного распределения после рассеяния, которые имеют те же моменты случайной величины, что и те, которые выдает теория Льюиса (Lewis). Проще говоря, теория Урбана делает попытку упрощенного описания теории Льюиса, которая основана на строгом решении уравнения диффузии. Подробно теорию Льиса можно изучить в его статье здесь.
2. G4WentzelVIModel - модель, относящаяся в смешанным алгоритмам. Такой алгоритм моделирует лобовые столкновения по отдельности используя модель однократного кулоновского рассеяния. Мягкие же столкновения с отклонением на небольшие углы, меньшие (по умолчанию, параметр можно менять) 0,2 рад, моделируются с помощью модели Wentzel-VI (как следует из названия).

Ионизационные потери протонов при прохождении через вещество. Потери энергии при прохождении через вещество протонов обусловлены прежде всего взаимодействием с электронными оболочками атомов, ионизацией и возбуждением их. Моделирование ионизационных потерь частицы делится на две фазы: дискретное моделирование каждого взаимодействия с рождением вторичных частиц и непрерывные потери энергии ниже задаваемого порога энергии частицы. Из имеющихся моделей в Geant4 для этого процесса нам подходит G4BetheBlochModel, название которой говорит само за себя.

• G4HadronElasticPhysics 

Этот раздел физики моделирует упругое столкновение протонов с ядрами вещества. Внутри него содержится нужная и единственная мадель для протонов по умолчанию G4ChipsElasticModel. Она моделирует упрогое рассеяние адронов на ядре, используя таблицы сечений и модель взаимодействия Chiral Invariant Phase Space (ChIPS). Подробнее о ней можно почитать в Geant4 Physics Reference Manual на странице 431. Вообще CHIPS моделирует фрагментацию адронных систем на адроны, кластеризацию нуклонов в ядре и механизмы взаимодействия адронов и фотонов с ядерной материей (как частный случай - упругое рассеяние). 

• G4HadronPhysicsQGSP_BIC ( и другие варианты: G4HadronPhysicsFTFP_BERT G4HadronPhysicsFTF_BIC и остальные классы, содержащие в названии BIC и BERT) или G4HadronPhysicsINCLXX

Раздел физики, отвечающий за моделирование неупругого взаимодействия адронов в ядрами. В нем существует несколько моделей и комбинаций моделей, однако, их замена не должна сильно влиять на результат эксперимента, так как неупругое столкновение протона с ядром - довольно редкое взаимодействие.

Прежде всего, первая часть аббревиатур QGS, FTF - обозначают то, что в соответствующих разделах физики используются модели Quark Gluon String и Fritiof, которые применимы только при энергиях протонов больше и значительно больше 1 ГэВ. Это не наш случай, поэтому формально эти модели в нашем эксперименте не будут вызываться вообще, поэтому первый кусок оббревиатуры не для нас, а для собратьев космофизиков.

Вторая часть аббревиатуры важнее. BIC - это Binary Cascade Model. В ней моделируется внутриядерный "путь" налетевшего протона и вторичных частиц такого неупругого взаимодействия. В каждый момент времени происходит взаимодействие налетевшего протона и только одного нуклона ядра, последовательно. Отсюда взялось название - бинарный каскад. Такой каскад останавливается как только энергия вторичных частиц становится ниже пороговой (она является нижней энергетической границей модели), которая составляет около сотни МэВ и может быть изменена. Дальше в игру вступают так называемая precompound model - модель, которая позволяет плавно перейти от кинетической стадии реакции к стадии равновесного ядра.

BERT - второй вариант окончания аббревиатуры - Bertini Cascade Model.  В ней ядро моделируется как слоистая концентрическая структура, каждый слой которой имеет постоянную плотность - имитация изменяющейся ядерной плотности. Каскад начинается, когда налетающая частица сталкивается с ядром и рождает вторичные частицы (при этом, заметьте, вторичнае частицы на этот момент являются кластерами нуклонов внутри ядра). Вторичные частицы в свою очередь могут взаимодействовать с другими нуклонами или отдать свою энергию ядру и быть поглощенными. Наступает фаза возбужденного ядра (описывается экситонной моделью), которое проходит стадию релаксации и испускает энергию в виде вторичных частиц. Эту фазу как раз и описывают G4PreCompoundModel и G4CascadeDeexcitation - две модели на выбор, которые используются внутри Bertini Cascade для моделирования переходя ядра в равновесное состояние. В Bertini Cascade по умолчанию стоит вторая, в Binary Cascade жестко поставлена первая. Это связано с особенностями описания модели. 

Кстати, подробнее о внутриядеоном каскаде и экситонной модели можно почитать здесь.