Занятие 2. Базовое описание материалов, геометрии и однородных полей в Geant4

Внимание! Данная статья подразумевает знание следующих понятий из C++:

1. Класс. Члены класса. Public, private члены класса.
2. Конструктор и деструктор класса.
3. Наследование класса и переопределение в дочернем функций материнского класса.
4. Объявление функций. Аргументы функций и значения, возвращаемые функциями.
5. Файл заголовка (header file) и файл с кодом (source file). Директивы  #include, #ifdef, #ifndef, #else, #endif.
6. Создание (инициализация) объекта определенного класса. Ключевое слово new и delete. Указатель на объект.
7. Ключевое слово static.

На примере Example1, который моделирует реальный физический эксперимент, мы изучим как задается геометрия в Geant4. Эксперимент схематично нарисован ниже.

Начнем с того, что изучим файл, содержащий функцию main, который называется Polyethylene.cc. В Geant4 подразумевается определенная структура программы, в которой происходит инициализация объектов в определенной последовательности. Если мы забежим вперед и, скажем, запустим расчет до того, как мы создали геометрию то, само собой, программа вылетит с ошибкой. Если говорить точнее, она даже не соберется.

В начале файла содержится большое количество необходимых инклюдов классов, которые нам будут нужны в дальнейшем. Пока стоит обратить внимание на следующий кусок:

#ifdef G4MULTITHREADED
#include "G4MTRunManager.hh"
#else
#include "G4RunManager.hh"
#endif

И ниже:

#ifdef G4MULTITHREADED
    G4MTRunManager* runManager = new G4MTRunManager;
    runManager->SetNumberOfThreads(1);
#else
    G4RunManager* runManager = new G4RunManager;
#endif

В первом куске содержится директива препроцессору, сообщающая о том, что если мы включили параметр G4MULTITHREADED на стадии cmake, то будет включен заголовок G4MTRunManager.hh (многопоточная версия класса G4RunManager), в другом случае будет включен G4RunManager. G4RunManager - самый главный класс в Geant4, который управляет всем: инициализацией геометрии, подключением физики, продвижением частицы через геометрию и т.д. Его функция Initialize() вызывается в main последней, уже после того как все условия эксперимента заданы.
Во втором куске инициализируется класс G4RunManager или, в многопоточной версии, класс G4MTRunManager с последующим заданием количества потоков.

Теперь мы дошли до этапа инициализации класса, содержащего геометрию.

    PolyethyleneDetectorConstruction* massWorld = new PolyethyleneDetectorConstruction;
    massWorld->RegisterParallelWorld(new PolyethyleneParallelWorld("PolyethyleneParallelWorld"));
    runManager->SetUserInitialization(massWorld);

В первой строчке объект класса с геометрией, собственно, инициализируется. Вторую строчку мы пока пропускаем - это предмет отдельного занятия. В третьей строчке указатель на объект класса с геометрией передается в G4RunManager. Теперь он знает все о нашей геометрии.
Теперь поговорим о том, что же это за класс PolyethyleneDetectorConstruction. Вы, наверно, обратили внимание, что его название явно связано с названием эксперимента. Это не случайно. Класс геометрии - это класс, который мы обязательно должны написать сами, наследуя класс G4VUserDetectorConstruction, содержащий пустую функцию Construct().
Давайте убедимся в этом сами. Наводим курсор на подсвеченное название класса PolyethyleneDetectorConstruction и нажимаем F2. Мы видим файл заголовка этого класса и вот такую строчку:

class PolyethyleneDetectorConstruction : public G4VUserDetectorConstruction

Эта конструкция и означает, что мы унаследовали класс G4VUserDetectorConstruction. Теперь если мы проведем ту же операцию с названием класса G4VUserDetectorConstruction, то окажемся в его файле заголовка и увидим такую строчку:

virtual G4VPhysicalVolume* Construct() = 0;

Она означает, что материнский класс G4VUserDetectorConstruction содержит пустую функцию Construct(), "= 0" всегда говорит о том, что функцию надо определить в дочернем классе обязательно.

Настало время задать материалы, которые вместе с геометрией описываются в функции Construct() нашего класса PolyethyleneDetectorConstruction. Для того, чтобы выделить подзадачу задания материалов отдельно (что является хорошим стилем в объектно ориентированном программировании) я создала свою собственную фунцию InitializeMaterials() в нашем классе в разделе private (так как никакому другому классу, кроме конкретно этого, эта функция не понадобится). Проверьте это сами. Обратите внимание, что функция ничего не возвращает, void.

Поговорим о том, как задаются материалы в Geant4. Теоретически, можно создать любые элементы, изотопы и материалы самим, вручную. Но обычно это не нужно, так как Geant4 содержит большую базу данных элементов и материалов. Класс, который обращается к этой базе данных (которая называется NIST) и передает указатели на нужные нам элементы и материалы называется G4NistManager. Это класс static и singleton, что означает что он доступен отовсюду из нашей программы и может быть инициализирован только один раз. Что означает, что может быть только один объект этого класса в нашей программе, и понятно почему - материал не может быть определен в двух разных местах, неоднозначно.
Итак, если класс G4NistManager static, то это означает, что мы можем прямо сейчас, находясь в функции InitializeMaterials() получить на него указатель. Смотрим в код файла PolyethyleneDetectorConstruction.cc (Напоминаю, что список всех файлов, относящихся к нашей программе, содержится слева в окне QtCreator).

void PolyethyleneDetectorConstruction::InitializeMaterials()
{
    G4NistManager* nistManager = G4NistManager::Instance();

Всегда, когда мы будем иметь дело со static файлом, мы будем получать на него указатель такоей же конструкцией через Instance() вместо new. new создает новый объект класса, а, повторюсь, объект этого класса должен быть только один. Теперь посомотрим на следующие строчки.

G4Element* H = nistManager->FindOrBuildElement(1);

Здесь мы получаем от G4NistManager указатель на необходимый элемент (в данном случае водород) путем вызова его функции FindOrBuildElement, которая принимает как аргумент номер элемента. Дальше идут точно такие же строчки, которые возвращают нам указатели на другие элементы, которые понадобятся нам при создании материалов. Обратите внимание, что мы не создаем элементы - мы всего лишь получаем на них указатели. Сами элементы уже созданы внутри G4NistManager, второй раз создавать их не имеет смысла.
После того, как мы получили указатели на все необходимые нам элементы, мы можем начать конструировать материалы. Опять обращу внимание, если мы используем материал из базы данных, то нам не надо его создавать через new, в противном случае мы будем использовать new.
Итак, первый пример получения указателя на материал из базы данных полностью аналогичен получению указателя на элемент из базы данных:

G4Material* Air = nistManager->FindOrBuildMaterial("G4_AIR");

Дальше я добавляю его в список материалов, который является private членом нашего класса PolyethyleneDetectorConstruction для того, чтобы мы могли воспользоваться этими материалами из другой функции нашего класса, а именно, из функции Construct(). Посмотрим еще ниже по коду и увидим такой кусок:

    G4Material* Brass = new G4Material("Brass", 8.50*g/cm3, 4);
    Brass->AddElement(Zn, 0.354);
    Brass->AddElement(Cu, 0.6175);
    Brass->AddElement(Pb, 0.025);
    Brass->AddElement(Fe, 0.0035);
    MaterialMap["Brass"] = Brass;

Здесь мы создаем материал латунь своими руками. Сделать так мне пришлось, потому что мне нужен был определенный сплав с определенной массовой долей каждого элемента и плотностью. Первой строчкой мы создаем новый объект класса G4Material через new. Перейдите по F2 в файл заголовка этого класса и просмотрите его функции в разделе public. Повторюсь, что к функциям из раздела private мы не имеем доступа извне. Можно увидеть что функция AddElement, которую мы вызываем впоследствии, принимает как аргумент указатель на элемент и его массовую долю в материале. В файле заголовка содержаться и другие варианты этой функции. Компилятор сам определяет какая функция имеется в виду по аргументам, которые вы в нее передаете.

Теперь перейдем к созданию геометрии и возвращаемся в функцию Construct(). Мы видим, что первой строчкой вызывается функция InitializeMaterials(). Пропускаем следующие две строчки, они касаются визуализации, и о них мы поговорим отдельно.
Теперь мы находимся на участке кода, который создает мировой объем. Это первое и обязательное, что мы дожны сделать при создании геометрии. Важно, чтобы все другие объемы, которые мы поместим внутрь мирового не выходили за его пределы, иначе программа выдаст ошибку.

    G4Box* world = new G4Box("World", 3*m, 3*m, 3*m);
    G4LogicalVolume *worldLogic = new G4LogicalVolume(world, MaterialMap["Air"], "WorldLogic");
    G4VPhysicalVolume *worldPhys = new G4PVPlacement(0, G4ThreeVector(), worldLogic, "WorldPhys", 0, false, 0);

Первой строчкой мы создаем solid. Это просто форма, которая содержит информацию только о собственном имени "World" и своих линейных размерах. Обратите внимание, что линейные размеры в Geant4 - это полуширины и радиусы. То есть, размер объекта отсчитывается от его центра. Итак, этой строчкой мы создали объект world класса G4Box с именем "World" и размерами 6х6х6 метров. С единицами измерения в Geant4 все просто, они пишутся через звездочку после цифры. Если в начале файла не стоит using namespace CLHEP, то вместо *m мы будем писать *CLHEP::m. CLHEP - файл, а если быть точным namespace, который содаржит определения и взаимосвязи между всевозможными единицами измерения.

Вторая строчка создает логический объем, который содержит информацию о материале, которым он заполнен, своем названии, и главное, содержит указатель на solid. Обратите внимание, что в конструктор класса G4LogicalVolume мы передаем указатель на solid - world. Получается, как будто, вложенная структура. И сейчас мы добавим еще один, третий слой.

G4VPhysicalVolume - класс, который сожержит информацию о расположении логического объема в пространстве, его повороте, материнском логическом объеме (в случае мирового объема - материнского объема нет), своем названии, и еще нескольких других параметрах. Каких - узнайте сами, нажав F2 и перейдя в файл заголовка. И главное, в конструкторе передается указатель на соответствующий ему логический объем. Обратите внимание, что конструктор G4PVPlacement отличасется от названия класса G4VPhysicalVolume. Это связано с тем, что класс G4PVPlacement - наследник класса G4VPhysicalVolume. Такова особенность языка C++. Для понимания этого нюанса вам будет необходимо продвинуться чуть дальше в понимании принципов объектно ориентированных языков и их качества - полиморфизма. Пока давайте просто запомним, что когда мы создаем физический объем, то используем конструктор G4PVPlacement.

Настало время обратиться к рисунку нашего эксперимента и обратить внимание, что весь объем состоит из повторяющихся структур, а именно тонкие листы латуни чередуются с толстыми слоями полиэтилена. Один лист латуни, заключенный между полуслоями полиэтилена составляет один канал. Поясню физический смысл канала. В данном эксперименте на слоистую структуру падал протонный пучок, и заряд, вследствие остановившихся протонов, накапливался в каналах. Лист латуни является проводником и был подключен к амперметру, который регистрировал ток. Такая конструкция наблюдалась в каждом канале. В итоге, к концу эксперимента у нас имелись данные о распределении заряда по глубине.
Начнем с того, что "нарисуем" один лист латуни с полуразмерами 7,5 х 7,5 х 0,00127 см. Вот как выглядел бы этот кусок кода, если бы нам нужно было нарисовать только один лист латуни:

    G4Box* brassSheet = new G4Box("BrassSheet", 7.5*cm, 7.5*cm, 0.00127*cm);
    G4LogicalVolume* brassLogic = new G4LogicalVolume(brassSheet, MaterialMap["Brass"], "BrassLogic");
    G4VPhysicalVolume* brassPhys  = new G4PVPlacement(0, G4ThreeVector(0,0,0), brassLogic, "BrassPhys", worldLogic, false, 0);

Обратите внимание, что вся конструкция полностью аналогична той, с помощью которой мы создавали мировой объем. Единственное отличием является то, что в строчке, где мы создаем физический объем для листа латуни, мы указываем материнский логический объем worldLogic. И еще, G4ThreeVector(0,0,0) и G4ThreeVector() значат одно и то же, вызывают конструктор вектора и создают нулевой вектор. Если мы используем нулевой вектор, то объем создастся в центре нашего мирового объема.
Первое, чем отличается код, который есть в программе от примера, представленного выше - это отсуствием определения типа объекта G4LogicalVolume* во второй строчке. Это связано с тем, что brassLogic является членом класса PolyethyleneDetectorConstruction и уже определен в его заголовке. Принято выделять переменные класса, называя их особенно. Я, например, называю их с большой буквы BrassLogic.
Второе отличие - отсутствие третьей строчки, которая создает физический объем. Нам необходимо сделать сборку повторяющейся структуры (один лист латуни, один лист полиэтилена) и отштамповать ее 65 раз, каждый раз сдвигаясь на ширину одного канала. Делаем мы это с помощью класса G4AssemblyVolume.

    G4AssemblyVolume* channelLogic = new G4AssemblyVolume();
    channelLogic->AddPlacedVolume(BrassLogic, firstPos, 0);
    channelLogic->AddPlacedVolume(CH2Logic, secondPos, 0);

Первая строчка - уже знакомая нам конструкция. Мы инициализируем объект channelLogic класса G4AssemblyVolume. Далее мы создаем "штамп" - добавляем два логических объема, меди и полиэтилена. firstPos и secondPos - это два вектора, отражающих расположение этих двух листов относительно друг друга.
Настало время отштамповать эту структуру 65 раз. Для этого мы создаем цикл:

    for( unsigned int i = 0; i <= 64; i++ )
    {
        G4ThreeVector offset(0, 0, (channelWidth*i));
        channelLogic->MakeImprint(worldLogic, offset, 0);
    }

Внутри цикла первая строчка рассчитывает смещение, вторая - вызвает функцию класса G4AssemblyVolume MakeImprint. Какие аргументы она принимает вы уже умеете выяснять сами.
Предпоследняя строчка в функции Construct() "рисует" последний лист латуни, который выбивается из повторяющегося рисунка. Мы используем только правую сторону из конструкции инициализации физического объема, так как указатель на этот физический объем нам нигде не понадобится. Мы вполне можем написать и целую строчку.

И, наконец, очень важная последняя строчка:

return worldPhys;

Она возвращает указатель на физический мировой объем, как и требует объявление функции Construct() в заголовке класса G4VUserDetectorConstruction:

virtual G4VPhysicalVolume* Construct() = 0;

В заголовке класса PolyethyleneDetectorConstruction:

G4VPhysicalVolume* Construct();

И в файле с кодом класса PolyethyleneDetectorConstruction:

G4VPhysicalVolume* PolyethyleneDetectorConstruction::Construct()

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

Теперь давайте поговорим о том, как задаются поля. Поля можно создавать какими угодно - однородными или изменяющимися по какому-либо закону, элетрические, магнитные, электро-магнитные, даже гравитационные или сильные. Все зависит от того, что мы заложим в класс поля G4Field, наследуя его.
Но сегодня, мы не будем ударяться в такие сложности (вы сможете в этом разобраться сами, когда понадобится, после того что мы сейчас разберем) и рассмотрим простой случай однородного магнитного поля с использованием стандартного класса G4UniformMagField, который является наследником класса G4Field.
В новом Geant4 10.01 есть очень важное правило. Все описания полей, а так же чувствительных детекторов (до которых мы еще доберемся) выделяются в отдельную функцию ConstructSDandField(). Эта функция тоже пустая в материнском классе G4VUserDetectorConstruction, но не содержит "= 0". Это означает, что ее определять не обязательно. То есть, если нам не нужны магнитные поля или чувствительные детекторы, мы вообще ее не объявляем в нашем классе геометрии.
Сделано это вот зачем. В новой версии Geant4  добавлена возможность многопоточного расчета. Грубо говоря, это возможность раскидать расчет на несколько параллельных вычислений с одной и той же геометрией и условиями. Объект геометрии остается общим для всех потоков. А вот поля и чувствительные детекторы "копируются" для каждого потока, и для каждого потока есть свой собственный объект поля и чувствительного детектора. Это нужно затем, что поля и детекторы активно участвуют в процессе трекинга частицы через геометрию, то есть в рассчете ее траектории. Поэтому для каждого потока, для каждой подзадачи, требуется свой собственный класс поля и детектора.
Функция ConstructSDandField() вызывается каждым потоком отдельно и создает для каждого потока его собственные классы поля и детектора. Если мы не используем многопоточность, то эта функция просто создает класс поля и детектора для нашего эксперимента.
Давайте посмотрим, что содержится в этой функции:

    MagneticField = new G4UniformMagField(G4ThreeVector(100500*gauss,100500*gauss, 0));
    G4FieldManager* fieldMgr = new G4FieldManager(MagneticField);
    BrassLogic->SetFieldManager(fieldMgr, false);
    CH2Logic->SetFieldManager(fieldMgr, false);
    //G4TransportationManager::GetTransportationManager()->SetFieldManager(fieldMgr);

В первой строчке мы создаем объект MagneticField класса G4UniformMagneticField. Обратите внимание, что объект я назвала с большой буквы, так как он является членом класса PolyethyleneDetectorConstruction, и его объявление содержится в заголовке этого класса.
Второй строчкой мы создаем объект fieldMgr класса G4FieldManager. Это класс, который должен содержать указатель на поле, используемое нами. При трекинге частицы, Geant4 обращается к нему.
Теперь у нас есть два варианта. Если мы хотим, чтобы магнитное поле у нас присутствовало только в каком-то логичиском объеме, а не во всем мире, то указатель на G4FieldManager, в который мы передали указатель на наше магнитное поле, мы передаем только  в нужный логический объем. Тогда этот FieldManager не является главным, а действует только в внутри объема, который соджержит на него указатель. (В программе может быть несколько объектов класса G4FieldManager). Попробуйте разные величины магнитного поля и запустите программу. Вы увидите, что треки протонов загибаютя по разному. А при значении магнитного поля больше 100000 Гаусс, протоны начинают вылетать из объема сбоку.
Теперь, если мы хотим задать магнитное поле во всем мировом объеме, то объект класса G4FieldManager, содержащий указатель на это магнитное поле нужно сделать главным. То есть передать в G4TransportationManager - static объект, который управляет прохождением частицы через геометрию. Этот объект единственный во всей программе и указатель на него мы получаем уже знакомым нам способом. Раскомментьте эту строчку и закоментьте две предыдущие. Запустите программу. Теперь вы видите, что магнитное поле действует во всем объеме.

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