Содержание |
Введение в использование GEANT4
Эта статья посвящена рассмотрению простейшего примера использования GEANT4, в котором происходит моделирование процесса рождения тормозного излучения. На рисунке приведена схема моделируемой установки (размеры даны в миллиметрах).
Установка состоит из вольфрамовой тормозной мишени толщиной 1мм и кремниевого детектора толщиной также 1 мм. Воздух из рассматриваемой области откачан. Вдоль оси Z на мишень подается пучок электронов с энергией 50 МэВ, который рождает в ней поток тормозного излучения. Задача состоит в измерении спектра фотонов с помощью детектора.
Проект на GEANT4, решающий данную задачу можно скачать по этой ссылке (http://download.neooffice.org/). Он состоит из следующих файлов:
ExampleG4.cc makefile vis.mac include src ./include: DetectorConstruction.hh PhysicsList.hh PrimaryGeneratorAction.hh SensitiveDetector.hh SteppingVerbose.hh ./src: DetectorConstruction.cc PhysicsList.cc PrimaryGeneratorAction.cc SensitiveDetector.cc SteppingVerbose.cc
Далее следует описание содержимого каждого файла.
Файл ExampleG4.cc
Программа на GEANT4 должна содержать определение нескольких основных классов, которые заключают в себя всю специфику конкретной задачи, и регистрацию этих классов в специальном объекте G4RunManager, который и управляет процессом моделирования. В число этих классов входят как обязательные:
- G4VUserDetectorConstruction, содержащий определение геометрии установки, и, обычно, определение используемых материалов и назначение чувствительных областей;
- G4VPhysicsList, подключающий моделирование интересующих физических процессов;
- G4VUserPrimaryGeneratorAction, описывающий источник первичных частиц в моделировании;
так и необязательные классы, без которых моделирование возможно: G4UserRunAction, G4UserEventAction и G4UserSteppingAction, позволяющие модифицировать поведение GEANT4 на том или ином этапе моделирования.
Определение каждого класса помещается в соответствующий отдельный файл, а в главном файле проекта ExampleG4.cc происходит сведение всего воедино и регистрация классов в G4RunManager.
Файл начинается с подключений заголовочных файлов. В Си++ объявления классов (то есть описание их полей, функций и констант) обычно размещаются в специальных заголовочных файлах с расширением .h или .hh, которые затем подключаются в файлы .cc директивой #include.
#include<G4RunManager.hh> #include<G4UImanager.hh> #include<G4UIterminal.hh> #include<G4VisExecutive.hh> #include<G4Material.hh> #include<G4UserRunAction.hh> #include<G4Run.hh> #include<iostream> #include<string> #include<CLHEP/Random/Random.h> #include<unistd.h> #include<time.h>
В списке подключаемых файлов сначала перечислены системные, входящие в состав GEANT4 и компилятора, а затем заголовочные файлы текущего проекта. Они отличаются способом задания имени: не в угловых скобках <>, а в кавычках. Эти файлы должны быть размещены в папке include.
#include "DetectorConstruction.hh" #include "PrimaryGeneratorAction.hh" #include "EventAction.hh" #include "SteppingAction.hh" #include "SteppingVerbose.hh" #include "PhysicsList.hh" using namespace std; const char macros[]="vis.mac";
Далее следует определение класса RunAction, который наследуется от класса G4UserRunAction и содержит функцию-член BeginOfRunAction, которая автоматически вызывается в начале каждого запуска. Здесь она просто выводит на экран порядковый номер запуска (0, 1 и т. д.). Для простоты определение этого класса не вынесено в отдельные файлы, а целиком приведено в ExampleG4.cc.
class RunAction: public G4UserRunAction
{
public:
void BeginOfRunAction(const G4Run* aRun)
{
G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
}
};
Далее следует определение функции main(). В Си++ main() является основной функцией программы, с которой начинается ее выполнение. Здесь происходит инициализация GEANT4 и все остальные вспомогательные действия, о которых говорилось выше.
int main(int argc,char** argv)
{
Установка экземпляра класса SteppingVerbose, который отвечает за печать подробной информации о каждом шаге частиц в процессе моделирования. Степень детализации выводимой информации зависит от числового параметра /stepping/verbose в файле vis.mac, 0 соответствует минимально подробной, а 9 максимально подробной информации о каждом шаге. Класс SteppingVerbose определен в файлах SteppingVerbose.hh и SteppingVerbose.cc. Оператор new создает объект — экземпляр этого класса и возвращает указатель на созданный объект.
G4VSteppingVerbose::SetInstance(new SteppingVerbose);
Настройка генератора случайных чисел. По умолчанию он возвращает одну и ту же последовательность случайных чисел, что удобно при отладке. Для генерирования более случайных последовательностей требуется задавать так называемое зерно (seed), которое в данном случае задается как сумма текущего времени в секундах и программного идентификатора. Это дает достаточно удовлетворительную случайность.
HepRandom::setTheSeed(time(0)+getpid());
Далее создается объект G4RunManager, который управляет запуском и остановкой моделирования.
G4RunManager * runManager = new G4RunManager;
Создается объект DetectorConstruction (см. класс DetectorConstrucion в файле DetectorConstruction.hh) и регистрируется в G4RunManager.
DetectorConstruction* detector_c = new DetectorConstruction;
runManager->SetUserInitialization(detector_c);
Так же создается и регистрируется в G4RunManager пакет физических процессов PhysicsList (см. соответствующие файлы).
G4VUserPhysicsList *p = new PhysicsList;
runManager->SetUserInitialization(p);
Создается и инициализируется объект класса G4VisExecutive, который позволяет различными способами визуализировать моделирование.
G4VisManager* visManager = new G4VisExecutive;
visManager->Initialize();
runManager->SetUserAction(new PrimaryGeneratorAction);
runManager->SetUserAction(new RunAction);
runManager->SetUserAction(new EventAction);
runManager->SetUserAction(new SteppingAction);
В конце концов вызывается метод G4RunManager::Initialize() и процесс инициализации GEANT4 завершается.
runManager->Initialize();
Печать информации о зарегистрированных материалах.
cout<<"===============================================================";
cout<<endl;
cout<< *(G4Material::GetMaterialTable()) << endl;
cout<<"===============================================================";
cout<<endl;
Наконец, через объект класса G4UImanager производится выполнение макрокоманд из файла vis.mac (на него указывает определенная выше переменная macros). Это удобно, потому что при изменении vis.mac не нужно перекомпилировать всю программу. В vis.mac находятся команды, непосредственно запускающие моделирование.
G4UImanager * UI = G4UImanager::GetUIpointer();
G4UIsession * session = new G4UIterminal();
UI->ExecuteMacroFile(macros);
Моделирование закончено. Освобождается память и программа завершается.
delete session;
delete visManager;
delete runManager;
return 0;
}
Класс DetectorConstruction
Геометрические характеристики моделирования в GEANT4 задаются через класс G4VUserDetectorConstruction. Чтобы сделать это, необходимо определить класс, наследующий G4VUserDetectorConstruction и определить в нем функцию Construct(), которая будет автоматически вызвана на этапе инициализации.
Объявление класса DetectorConstruction находится в файле DetectorConstruction.hh.
DetectorConstruction.hh
В начале файла находится так называемый include guard, необходимый во всех заголовочных файлах.
#ifndef DetectorConstruction_h #define DetectorConstruction_h 1
Подключение заголовочных файлов GEANT4. В GEANT4 используется следующее соглашение: каждый класс объявляется в отдельном файле, и имя файла совпадает с именем класса. Файл globals.hh содержит физические константы и т. п.
#include<globals.hh> #include<G4VUserDetectorConstruction.hh> #include<G4VSolid.hh> #include<G4LogicalVolume.hh> #include<G4VPhysicalVolume.hh> #include<G4Material.hh>
Для удобства определен класс World, который соответствует внешнему объему в геометрии GEANT4, внутри которого находится вся моделируемая установка. Этот объем можно было бы создать и просто через цепочку Solid→Logic→Physic, но можно и инкапсулировать ее в одном классе World.
class World {
protected:
G4VSolid *solid;
G4LogicalVolume *logic;
G4VPhysicalVolume *physic;
G4Material *mater;
double sizex, sizey, sizez;
public:
World(double size_x, double size_y, double size_z, G4Material *mater_=NULL);
operator G4LogicalVolume*() {return logic;}
G4LogicalVolume *getLogic() {return logic;}
// void setLogic(G4LogicalVolume *volA);
G4VSolid *getSolid() {return solid;}
G4VPhysicalVolume *getPhysic() {return physic;}
};
Объявление класса DetectorConstruction, наследуемого от G4VUserDetectorConstruction.
class DetectorConstruction : public G4VUserDetectorConstruction
{
Открытые члены класса DetectorConstruction. В Си++ бывает три уровня доступа к членам класса:
- public, когда любой код может обратиться к данному члену (изменить
или прочитать значение переменной или вызвать функцию-член);
- private: когда только функции данного класса могуть обращаться
к данному члену;
- protected: когда к данному члену могут обращаться функции
данного класса и классов, наследованных от данного.
Далее идет объявление конструктора и деструктора класса DetectorConstruction — функций, которые автоматически вызываются при создании и уничтожении объекта. В частности, если конструктор будет объявлен закрытым (private или protected), то создать объект будет невозможно.
public:
DetectorConstruction();
~DetectorConstruction();
Функция Construct заключает в себе основную функциональность класса DetectorConstruction. Она создает геометрию и материалы.
G4VPhysicalVolume* Construct();
protected:
Также объявляется защищенная переменная-указатель на World. Она будет проинициализированна позднее.
World *world;
В конце объявления каждого класса по правилам синтаксиса Си++ должна стоять точка с запятой.
}; #endif
Определение объявленного в файле DetectorConstruction.hh класса DetectorConstruction находится в файле DetectorConstruction.cc.
DetectorConstruction.cc
#include<G4NistManager.hh> #include<G4Box.hh> #include<G4Tubs.hh> #include<G4LogicalVolume.hh> #include<G4PVPlacement.hh> #include<G4SDManager.hh> #include<G4VisAttributes.hh> #include"DetectorConstruction.hh" #include"SensitiveDetector.hh" using namespace std;
Здесь определяется макрос Mat(), просто для того чтобы сократить запись: теперь вместо того чтобы писать G4NistManager::Instance()->FindOrBuildMaterial(«G4_Si») достаточно написать Mat("G4_Si"). При этом будет произведена инициализация соответствующего материала в базе данных GEANT4 (если он в ней содержится, конечно). Список имеющихся в ней материалов можно посмотреть в файле $G4INSTALL/source/material/management/src/G4NISTMaterialBuilder.cc. Если нужный материал не содержится в базе, то его можно определить самостоятельно, но рассмотрение этого вопроса выходит за рамки статьи.
#define Mat(x) (G4NistManager::Instance()->FindOrBuildMaterial(x))
Определение конструктора класса World. Он принимает четыре параметра: ширину, высоту, глубину и материал, и создает в качестве материнского объема бокс с этими параметрами через цепочку вызовов Solid→Logic→Physic, о которой немного позже.
World::World(double size_x, double size_y, double size_z, G4Material *mater_):
mater(mater_), sizex(size_x), sizey(size_y), sizez(size_z)
{
solid = new G4Box("world", sizex/2, sizey/2, sizez/2);
logic = new G4LogicalVolume( solid, mater, "World", 0, 0, 0);
physic = new G4PVPlacement(0, G4ThreeVector(), logic, "World",
0, false, 0);
}
Конструктор и деструктор DetectorConstruction пустые, т. к. и сам класс очень простой.
DetectorConstruction::DetectorConstruction()
{
}
DetectorConstruction::~DetectorConstruction()
{
}
Основную работу в DetectorConstruction выполняет функция Construct()
G4VPhysicalVolume* DetectorConstruction::Construct()
{
Для начала создается материнский объем размером 30x30x30 см., «заполненный» вакуумом.
world = new World(30*cm, 30*cm, 30*cm, Mat("G4_Galactic"));
Затем создается тормозная мишень. Каждый элемент геометрии GEANT4 задается тремя объектами. G4VSolid описывает его геометрические свойства: куб, сфера и т. п. В данном случае это бокс G4Box, то есть параллелепипед. Его параметры это половинные размеры.
G4Box *solidTgt = new G4Box("solidTgt", 2.5*cm, 2.5*cm, 0.5*mm);
На следующем этапе создается логический объем G4LogicalVolume, который содержит информацию о материале и магнитных свойствах среды. В данном случае никаких магнитных свойств нет.
G4LogicalVolume *logiclTgt = new G4LogicalVolume(solidTgt,
Mat("G4_W"),
"logiclTgt");
Третий этап — это физический объем, в котором содержится информация о положении объекта относительно материнского, то есть World. Здесь оно задается вектором G4ThreeVector(0,0,-5*cm), то есть -5 см по оси Z.
G4PVPlacement *physilTgt = new G4PVPlacement(0,
G4ThreeVector(0,0,-5*cm),
logiclTgt,
"physilTgt",
world->getLogic(),
false,
0);
Аналогично создается объем кремниевого детектора.
G4Box *solidDet = new G4Box("solidDet", 2.5*cm, 2.5*cm, 0.5*mm);
G4LogicalVolume *logicDet = new G4LogicalVolume(solidDet,
Mat("G4_Si"),
"logicDet");
G4PVPlacement *physiDet = new G4PVPlacement(0,
G4ThreeVector(0,0,5*cm),
logicDet,
"physiDet",
world->getLogic(),
false,
0);
Для того, чтобы детектор реагировал на пролет частиц, он должен быть назначен так называемой чувствительной областью. Для этой цели в G4LogicalVolume есть специальное поле SensitiveDetector, которое содержит указатель на объект класса G4VSensitiveDetector. Этот объект автоматически вызывается каждый раз, когда очередной шаг моделирования частицы попадает внутри данного объема. Соответствующим образом программируя класс SensitiveDetector можно получать выходные данные моделирования (см. раздел SensitiveDetector далее).
SensitiveDetector *detector = new SensitiveDetector("hi there");
Объект SensitiveDetector должен быть зарегистрирован в G4SDManager.
G4SDManager* SDman = G4SDManager::GetSDMpointer();
SDman->AddNewDetector(detector);
Здесь SensitiveDetector сопоставляется кремниевому детектору.
logicDet->SetSensitiveDetector(detector);
Стенки объекта World делаются прозрачными, чтобы не мешали на визуализации.
world->getLogic()->SetVisAttributes (G4VisAttributes::Invisible);
В конце концов успешно выполнившаяся функция Construct должна возвратить указатель на физический объем объекта World.
return world->getPhysic(); }
SenstiveDetector
Объекты SensitiveDetector обрабатывают информацию о каждом шаге моделирования внутри области геометрии, которой они назначены.
SensitiveDetector.hh
#ifndef SENSITIVEDETECTOR #define SENSITIVEDETECTOR #include<G4VSensitiveDetector.hh> class G4Step; class G4TouchableHistory;
Здесь объявляется класс SensitiveDetector. В данном случае он должен строить распределение выделенной энергии в детекторе. Рассчитывается гистограмма, а затем значения из ее столбцов записываются в текстовый файл, чтобы затем с ними работать в программе построения графиков вроде gnuplot, Excel или Origin.
class SensitiveDetector: public G4VSensitiveDetector
{
private:
Число столбцов в гистограмме.
static const int NOBINS = 1000;
Максимальный и минимальный пределы графика. Значения этих констант задаются в файле SensitiveDetector.cc.
const double HIST_MAX;
const double HIST_MIN;
Гистограмма будет представлена массивом int.
int histogram[NOBINS];
public:
SensitiveDetector(G4String name);
~SensitiveDetector();
G4bool ProcessHits(G4Step *step, G4TouchableHistory *hist);
void EndOfEvent(int nEvent);
};
#endif /* SENSITIVEDETECTOR */
SensitiveDetector.cc
В начале все как обычно: подключаются используемые библиотечные классы.
#include<G4Step.hh> #include<fstream> #include<iostream> #include"SensitiveDetector.hh" using namespace std;
Конструктор класса SensitiveDetector. Он принимает в качестве параметра имя, которое имеет тип G4String. При создании объекта это записывается так: new SensitiveDetector(«a name»). Кроме того здесь же инициализируются константы, задающие верхний и нижний предел на графике. Почему это нельзя было сделать в SensitiveDetector.hh? Можно, но для смеха я поместил этот код сюда. Потом обнуляются ячейки массива histogram, которые по умолчанию содержат мусор.
SensitiveDetector::SensitiveDetector(G4String name): G4VSensitiveDetector(name),
HIST_MAX(10*MeV),
HIST_MIN(0 *MeV)
{
for(int i = 0; i<NOBINS; i++)
histogram[i] = 0;
}
Основная функция этого класса — ProcessHits(). Всякий раз, когда очередной шаг моделирования попадает в объем, которому принадлежит данный SensitiveDetector, вызывается эта функция.
G4bool SensitiveDetector::ProcessHits(G4Step *step, G4TouchableHistory *hist)
{
Затем мы получаем полную энергию частицы.
double energy = step->GetTrack()->GetDynamicParticle()->GetTotalEnergy();
Теперь мы заносим полученное значение в гистограмму. Эта задача сводится к увеличению на 1 столбца гистограммы, соответствующего данной энергии.
double bin_width = (HIST_MAX - HIST_MIN) / NOBINS;
int index = int(floor((energy_deposit-HIST_MIN)/bin_width));
if(index >= 0 && index < NOBINS)
histogram[index]++;
В последнюю очередь мы уничтожаем данную частицу, чтобы не регистрировать ее дважды. При этом частица выводится из моделирования. Так как нам не интересна ее дальнейшая судьба, то такой простой подход уместен. В противном случае, если нужно моделировать слои детекторов и т. п., может понадобиться более сложная схема, когда происходит разделение частиц по TrackID и тому подобные вещи.
step->GetTrack()->SetTrackStatus(fStopAndKill);
return true;
}
Деструктор ~SensitiveDetector. Эта функция вызывается автоматически при удалении объекта, поэтому здесь удобно разместить вывод результатов в файл.
SensitiveDetector::~SensitiveDetector()
{
В Си++ работа с файлами происходит так же, как и с экраном, через потоки. Создаем поток, вывод которого направляется в файл spectrum.dat.
std::ofstream file("spectrum.dat");
И записываем в него гистограмму в формате «энергия—число отсчетов».
double bin_width = (HIST_MAX - HIST_MIN) / NOBINS;
for(int i = 0; i<NOBINS; i++)
{
double energy = i*bin_width + HIST_MIN;
file << std::setw(15) << energy/MeV << " "
<< std::setw(15) << histogram[i] << std::endl;
}
}
PhysicsList
В GEANT4 все используемые в моделировании физические процессы должны быть заранее подключены. Делается это в классе PhysicsList, то есть пакет физики. В нем должны создаваться определения всех частиц и для каждой частицы должны выбираться процессы, которые с ней могут происходить. Но так как создание такого пакета вручную довольно сложное дело, где надо учесть массу деталей, обычно используются готовые пакеты. В данном примере это так называемый стандартный PhysicsList, используемый в большинстве программ на GEANT4. Он включает в себя следующие процессы: для фотонов
- фотоэффект;
- комптоновское рассеяние;
- рождение пар.
Для электронов, позитронов и мюонов:
- тормозное излучение;
- множественное рассеяние;
- ионизация среды;
- аннигиляция с античастицами.
Кроме того для всех типов частиц подключается процесс Transportation, отвечающий за перемещение частиц в пространстве с учетом влияния магнитного поля.
Код класса PhysicsList находится в файлах PhysicsList.hh и PhysicsList.cc. Он совершенно стандартный, и единственный параметр, который в нем может потребоваться изменить — это так называемый кат. Он задается в PhysicsList.cc:
PhysicsList::PhysicsList(): G4VUserPhysicsList()
{
defaultCutValue = .1*cm;
SetVerboseLevel(1);
}
Как видно, это величина с размерностью длины. Смысл ее таков: для каждого типа частиц — электронов, фотонов и т. д.,— и для каждого материала вычисляется энергия, при которой эта величина совпадает со средним пробегом в данной среде. Затем, в ходе моделирования, если на некотором шаге должна родиться вторичная частица, но ее энергия меньше, чем энергия ката, то эта частица не рождается, а считается поглощенной в среде в данной точке. При этом ее энергия добавляется к TotalEnergyDeposit. Например, если в тормозной мишени из вольфрама фотон с энергией 10 МэВ рождает пару электрон-позитрон, а величина ката равна 1 мм, то эти частицы не будут рождаться, т. к. энергия каждой 5 МэВ и их средний пробег в вольфраме меньше 1 мм. Понятно, что
1. во-первых, кат должен быть меньше или равен характерной
толщине моделируемого объекта;
2. а во-вторых, что значение ката влияет прежде всего на низкоэнергетическую
часть спектра.
В примере стоит 1 мм, что соответствует толщине мишени и детектора.
PrimaryGeneratorAction
PrimaryGeneratorAction отвечает за создание первичной частицы, с которой начинается процесс моделирования события. Первичной частицей может быть электрон из ускорителя, фотон, излученный из радиоактивного источника, или вообще любая другая частица.
PrimaryGeneratorAction.hh
#ifndef PrimaryGeneratorAction_h #define PrimaryGeneratorAction_h 1 #include <G4VUserPrimaryGeneratorAction.hh> class G4ParticleGun; class G4Event;
Класс PrimaryGeneratorAction должен наследоваться от G4VUserPrimaryGeneratorAction.
class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
{
public:
PrimaryGeneratorAction();
~PrimaryGeneratorAction();
public:
void GeneratePrimaries(G4Event*);
private:
Задается объект G4ParticleGun, который позволяет ставить частицу в любой точке.
G4ParticleGun* particleGun; }; #endif
PrimaryGeneratorAction.сс
#include "PrimaryGeneratorAction.hh" #include<G4Event.hh> #include<G4ParticleGun.hh> #include<G4ParticleTable.hh> #include<G4ParticleDefinition.hh> #include<globals.hh>
Конструктор класса. Здесь удобно создать G4ParticleGun и настроить его на определенный тип частиц. В данном случае это электроны с энергией 50 МэВ, а точка старта -10 см по оси Z. Кроме того задается направление импульса вдоль Z, а сам импульс вычисляется автоматически.
PrimaryGeneratorAction::PrimaryGeneratorAction()
{
G4int n_particle = 1;
particleGun = new G4ParticleGun(n_particle);
G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
G4ParticleDefinition* particle = particleTable->FindParticle("e-");
particleGun->SetParticleDefinition(particle);
particleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.));
particleGun->SetParticlePosition(G4ThreeVector(0, 0, -10*cm));
particleGun->SetParticleEnergy(50*MeV);
}
В деструкторе G4ParticleGun удаляется.
PrimaryGeneratorAction::~PrimaryGeneratorAction()
{
delete particleGun;
}
GeneratePrimaries — это основная функция класса PrimaryGeneratorAction. Она вызывается в начале каждого события, чтобы сгенерировать первичные частицы, с которых начинается моделирование. Параметры частиц были раз и навсегда заданы в конструкторе, а здесь производится запуск ParticleGun.
void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
particleGun->GeneratePrimaryVertex(anEvent);
}
Компиляция и запуск
Компиляцией программы управляет мейкфайл. Здесь необходимо придерживаться соглашения GEANT4 о том, что имя программы, указанное на строчке name, должно совпадать с именем главного файла (ExampleG4.cc). Скомпилированная программа будет также в этом случае называться ExampleG4.
name := ExampleG4 G4TARGET := $(name) G4EXLIB := true G4DEBUG=1 CPPFLAGS+=-Wno-overloaded-virtual -Wno-unused .PHONY: all all: lib bin include $(G4INSTALL)/config/binmake.gmk
Для полной сборки проекта надо выполнить в корневой папке команду make clean && make. В дальнейшем, если модификации файлов были незначительными, можно не запускать make clean, а писать просто make.
Программа запускается командой $G4BINDIR/ExampleG4. После запуска вы увидите приблизительно такое окно, в котором будут отображаться траектории частиц.
