На мой взгляд, это самое сложное, с чем может столкнуться пользователь при моделировании своего эксперимента. Тем не менее, я постараюсь изложить данный этап моделирования как можно более понятно и по возможности подробно. Повторюсь, что если у читателя возникнут вопросы, я всегда готова ответить на них в комментариях.
Как и в предыдущем посте про физические модели, сегодня мы будем говорить о физике применительно к адронной терапии. Примерный набор физических процессов, которые используются в этой области физики можно посмотреть в примере "hadrontherapy" в директории "/examples/advanced/hadrontherapy" внутри дистрибутива Geant4 10.01 (файл HadrontherapyPhysicsList.сс). Или в нашем примере - там, по сути, используется та же физика.
Теперь поговорим подробно о способах, которыми эта физика задается. Всю систему моделирования физики внутри Geant4 можно разделить на три принципиальных элемента:
Как и в предыдущем посте про физические модели, сегодня мы будем говорить о физике применительно к адронной терапии. Примерный набор физических процессов, которые используются в этой области физики можно посмотреть в примере "hadrontherapy" в директории "/examples/advanced/hadrontherapy" внутри дистрибутива Geant4 10.01 (файл HadrontherapyPhysicsList.сс). Или в нашем примере - там, по сути, используется та же физика.
Теперь поговорим подробно о способах, которыми эта физика задается. Всю систему моделирования физики внутри Geant4 можно разделить на три принципиальных элемента:
- Физическая модель - генерирует конечное состояние частицы после взаимодействия и обновляет параметры моделируемой частицы.
- Физический процесс - содержит необходимые таблицы поперечных сечений и используемые физические модели.
- Список физических процессов (Physics List) - как понятно из названия, элемент, который содержит в себе набор физических процессов для каждой частицы.
При написании программы пользователь обязательно должен создать свой набор физических процессов в своем классе, который должен быть унаследован от G4VUserPhysicsList. После чего указатель на него он передает в G4RunManager. Вот как это должно выглядеть в коде в функции main:
Такие готовые наборы физики подключаются очень просто. Вот так:
Второй по убыванию сложности метод задания физики - это наследование класса G4VModularPhysicslist, который, в свою очередь является наследником класса G4VUserPhysicsList. У нас получается две ступени наследования!
В функции main это будет выглядеть так:
Итак, что у нас есть: класс 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()? Обратите внимание на следующий кусок кода:
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. И не используйте нулей, код программы игнорирует модели под нулевым номером. Странно, согласна.
В завершение я хочу упомянуть, что у электромагнитных процессов есть механизмы ограничения шага частицы, которые повышают точность моделирования. Они предусмотрены и для торможения частицы, и для многократного рассеяния, и общие. Вообще, физику можно использовать и без настройки этих параметров, так как они разумно выставлены по умолчанию. Однако, в специфических экспериментах это может понадобится. По поводу этих алгоритмов я напишу отдельный пост, но, к сожалению, на английском языке. Дело в том, что это может быть интересно иностранному сообществу пользователей тоже.