Geant4 무작정 따라하기 시리즈의 열한번째. Geant4에서 스코어링을 구현하는 방법에 대해 알아봅니다.
지난 글에서 살펴본 스코어링 이론을 바탕으로, 이번 글에서는 스코어링을 직접 구현해보도록 하겠습니다.
스코어링을 구현하는 방법은 크게 다음의 3단계를 통해 진행됩니다.
차례차례 살펴보도록 합시다.
Sensitive Detector를 설계하는 방법은 지난 글에서 말씀드렸습니다. 이 시리즈에서는 Multi Functional Detector(MFD)에 사용자가 원하는 Primitive Scorer들을 Register하여 설계하는 방법만 다루겠다고 하였지요.
이 내용을 코드로 입력하는 곳도 정해져 있습니다. 바로 DetectorConstruction.cc 파일에 있는 ConstructSDandField() 함수 내부입니다. DetectorConstruction.cc 파일을 열어보시면, 대부분의 경우 맨 아래쪽 부근에 다음과 같은 함수가 있을 것입니다. 이 함수의 중괄호({}) 내부에 작성할 것입니다.
물론, 관련한 헤더는 당연히 DetectorConstruction.cc 파일의 맨 위에 포함시켜 주셔야 합니다. 여기서 필요한 헤더는 다음과 같습니다.
Multi Functional Detector를 만드는 과정은 딱 두 가지만 하시면 됩니다.
코드로는 다음의 두 줄입니다.
1auto mfd = new G4MultiFunctionalDetector("Detector");
2G4SDManager::GetSDMpointer()->AddNewDetector(mfd);
mfd라는 변수명으로 객체를 만들고, 이를 추가해주었습니다.
좀 더 자세히 살펴보도록 하죠.
Geant4에서는 Multi Functional Detector를 만들기 위해 이름 그대로 G4MultiFunctionalDetector라는 클래스를 제공하고 있습니다. 사용자는 이 클래스의 생성자를 이용하여 객체를 한 개 만들어주면 됩니다.
이 클래스의 생성자는 다음과 같습니다.
1G4MultiFunctionalDetector(G4String name);
문자열 형태의 입력값 한 개만 자유롭게 넣어주면 됩니다. 입력 인자로 넣어준 이 값은 MFD의 이름이 되는데, 나중에 HCofThisEvent에서 사용되니 적절한 이름으로 잘 지어주세요. 위 예시 코드에서는 “Detector"라는 이름을 지어주었습니다.
Geant4에서 Sensitive Detector를 사용하기 위해서는 반드시 G4SDManager에 그 SD를 추가해주어야 합니다. 그래야 Geant4가 입자를 수송하면서 SD 목록 중 조건에 부합하는 것들에 대해 스코어링을 수행해주기 때문입니다.
G4SDManager는 Geant4가 제공하는 것을 가져다 쓰면 되며, 이 클래스가 제공하는 AddNewDetector()라는 함수를 통해 우리가 설계한 SD를 추가하게 됩니다. 이 함수의 원형은 다음과 같습니다.
1void AddNewDetector(G4VSensitiveDetector *aSD);
이 함수의 입력 인자로, 앞서 저희가 만들어둔 MFD 객체를 그대로 써주시면 됩니다.
이제 MFD는 만들었으니, 실제 스코어링 기능을 수행하는 Primitive Scorer를 등록해야 합니다. 이 내용도 ConstructSDandField() 함수 내부에 위에서 작성한 두 줄에 이어서 계속 써주시면 됩니다.
여기도 딱 두 가지만 해주시면 됩니다.
예를 들어 해당 지오메트리 내에 deposit된 에너지의 총합을 스코어링하는 G4PSEnergyDeposit을 등록하고 싶다면, 다음과 같이 두 줄을 작성하면 됩니다.
1auto psEDep = new G4PSEnergyDeposit("EDep");
2mfd->RegisterPrimitive(psEDep);
psEDep이라는 변수명으로 G4PSEnergyDeposit 클래스의 객체를 만들고, 이를 앞서 만들어둔 mfd 변수명으로 정의된 MFD 객체에 등록해주었습니다.
좀 더 자세히 살펴보겠습니다.
Primitive Scorer의 종류가 너무 많기 때문에, 이 글에서 모든 PS의 생성자를 살펴볼 수는 없습니다. 하지만 모든 PS의 생성자에 있어 공통적으로 중요한 내용이 있습니다. 대표적으로 G4PSEnergyDeposit 클래스의 생성자를 살펴봅시다.
1G4PSEnergyDeposit(G4String name, G4int depth=0);
2G4PSEnergyDeposit(G4String name, const G4String& unit, G4int depth=0);
MFD 객체를 만들 때와 마찬가지로, 문자열 형태의 입력값 한 개만 자유롭게 넣어주면 됩니다. 이 입력 인자는 어느 PS를 사용하든 무조건 입력해주도록 되어있으며, 이 값이 해당 PS의 이름이 됩니다.
이 또한 나중에 HCofThisEvent에서 사용되니 적절한 이름으로 잘 지어주세요. 위 예시 코드에서는 “EDep"라는 이름을 지어주었습니다.
MFD에 PS를 등록할 때에는, G4MultiFunctionalDetector 클래스가 제공하는 RegisterPrimitive() 함수를 이용합니다. 이 함수의 원형은 다음과 같습니다.
1G4bool RegisterPrimitive(G4VPrimitiveScorer *aPS);
이 함수의 입력 인자로, 앞서 저희가 만들어둔 PS의 객체를 그대로 써주시면 됩니다.
하나의 MFD에 여러 개의 PS를 등록하는 방법도 간단합니다. 위에서 설명한 두 과정만 계속 반복해주시면 됩니다.
예를 들어 MFD를 만들고 G4PSEnergyDeposit과 G4PSDoseDeposit 두 개의 PS를 등록하는 경우를 코드로 작성한다면 다음과 같습니다.
1auto mfd = new G4MultiFunctionalDetector("Detector");
2G4SDManager::GetSDMpointer()->AddNewDetector(mfd);
3auto psEDep = new G4PSEnergyDeposit("EDep");
4mfd->RegisterPrimitive(psEDep);
5auto psDoseDep = new G4PSDoseDeposit("DoseDep");
6mfd->RegisterPrimitive(psDoseDep);
Filter는 PS의 추가옵션 같은 기능입니다. Filter를 세팅하기 위해서는 PS를 등록할 때 Filter 관련 내용을 추가해야 합니다.
예를 들어, G4PSEnergyDeposit을 이용하는데, 전자(e-)가 deposit한 값만 스코어링하고 싶다면, 다음과 같이 작성하면 됩니다.
1auto psEDep = new G4PSEnergyDeposit("EDep");
2auto filterElectron = new G4SDParticleFilter("e-Filter", "e-"); // Create filter object
3psEDep->SetFilter(filterElectron); // Set the filter to PS
4mfd->RegisterPrimitive(psEDep);
각 Filter마다 생성자의 형태나 사용방법은 상이하므로, 여기서 자세히 다루지는 않겠습니다.
이제 관심 있는 지오메트리의 Logical Volume에, 앞서 설계한 SD를 세팅해줄 차례입니다. 이 내용 또한 ConstructSDandField() 함수 내부에 작성하며, 아주 간단하게 한 줄이면 끝납니다. SetSensitiveDetector() 함수를 사용하면 됩니다.
예를 들어, “phantom"이라는 이름을 가진 Logical Volume에 앞서 만든 mfd 변수명을 가진 SD를 세팅한다면, 다음과 같이 입력하면 됩니다.
1SetSensitiveDetector("phantom", mfd);
좀 더 자세히 살펴보도록 합시다. 이 함수의 원형은 다음과 같습니다.
1void SetSensitiveDetector(const G4String& logVolName,
2 G4VSensitiveDetector* aSD, G4bool multi = false);
3void SetSensitiveDetector(G4LogicalVolume* logVol, G4VSensitiveDetector* aSD);
SetSensitiveDetector 함수를 사용할 때, 관심 있는 지오메트리의 Logical Volume을 입력하는 방법이 두 가지입니다.
일단 이 시리즈에서 추천하는 바는 그냥 Logical Volume의 이름을 입력하는 것입니다. 제가 이전에 Logical Volume을 정의하는 방법에 대해 적은 글에서, Logical Volume의 이름을 지을 때 다른 Logical Volume과 겹치지 않게끔 고유의 이름을 권장한다고 하였습니다. 그 이유 중 하나가 이것입니다. 이름이 중복되면 그 중 어느 지오메트리에 SD를 세팅할 지가 불분명하기 때문입니다.
물론, 이름이 중복되는 모든 Logical Volume들에게 다 동일한 SD를 세팅하고자 한다면, 이 함수의 multi 인자 값을 true로 설정하여 사용할 수도 있습니다.
이 시리즈에서는 지역변수나 변수의 범위(scope) 등을 설명하기에는 지면이 부족하여, 그냥 Logical Volume의 이름을 사용하는 방법만 설명합니다.
일반적으로 Logical Volume 객체는 Construct() 함수 내에 지역변수로 만들기 때문에, ConstructSDandField() 함수에서는 그 객체를 호출하지 못하기 때문이지요.
물론 객체를 불러오기 위한 몇 가지 방법들이 있습니다만 여기서는 생략하도록 하겠습니다.
여기까지의 내용은 모두 DetectorConstruction.cc 파일에 있는 ConstructSDandField() 함수 내부에 작성하였습니다. 모두 마무리하셨다면, 일단 SD를 설계하고 이를 Logical Volume에 세팅하는 것까지 마친 상태입니다.
예를 들어 G4PSEnergyDeposit 클래스를 이용하여 설계한 SD를 “phantom"이라는 이름을 가진 Logical Volume에 세팅한다면, 다음과 같이 코드가 작성되어야 합니다.
여기서 중간 점검을 하는 이유는, 다음 내용부터는 작성하는 파일이 달라지기 때문입니다. 놓치지 말고 잘 따라오세요.
이제 마지막 단계입니다. 앞서 두 단계를 통해 Geant4는 알아서 매 Event마다 Hit을 생성하고 이를 HitsCollection에 담아 HCofThisEvent로 묶어서 가지고 있을 것입니다. Event가 끝날 때 사용자가 이 Hit을 꺼내서 확인하지 않으면, 그 정보들은 그대로 사라집니다. 메모리를 날려버리고 새로운 다음 Event에 대한 정보를 기록하지요.
사용자가 할 일은 Event가 끝날 때마다 Hit을 가져와서 확인하는 것입니다.
Geant4는 매 Event마다 특정 행동을 반복적으로 수행하는 용도로 활용할 수 있게끔 하기 위해, EventAction이라는 개념을 제공하고 있습니다. 이 시리즈에서는 EventAction.cc라는 파일을 열어서 이용하시면 됩니다. 여기에는 BeginOfEventAction()이라는 함수와 EndOfEventAction()이라는 함수가 있는데, 이름 그대로 매 Event가 시작되기 직전과 끝난 직후에 해당 함수가 실행되는 식으로 동작합니다. 즉, 우리처럼 매 Event가 끝날 때마다 무언가를 하고 싶다면 EndOfEventAction() 함수 안에 할 일을 적으면 되는 것이죠.
EventAction.cc 파일을 열어보시면 아래쪽에 EndOfEventAction() 함수가 있을 것입니다. 이 함수 안에 내용을 작성할 것입니다.
이미 다음과 같은 내용이 작성되어 있네요.
1auto HCE = anEvent->GetHCofThisEvent();
2if (!HCE)
3 return;
의미는 그냥 써있는 그대로 읽으시면 됩니다. 한글로 쓰자면 다음과 같습니다.
어때요, 참 쉽죠?
HCE 변수가 유효하지 않다면 if 조건문에 의해 함수가 종료되어 버리므로, 이 내용 아래 부분에 있는 내용은 반드시 HCE가 유효한 경우에만 실행될 것입니다.
이제 HCofThisEvent로부터 Hit까지 가보도록 합시다.
HCofThisEvent는 HitsCollection들의 묶음이라고 하였습니다. 이 때 Geant4는 수많은 HC들을 쉽게 구분하기 위해 HC마다 0 이상의 양의 정수값으로 된 고유의 ID번호를 부여합니다. 다만, 사용자가 이 정수값을 직접 기억할 필요는 없습니다. 사용자는 HC의 이름을 이용하여 이 ID번호를 찾아올 수 있기 때문입니다.
여기서 HC의 이름은 Sensitive Detector를 설계할 때 결정됩니다. 우리가 했던 것처럼 MFD와 PS를 이용할 경우에는, HC의 이름이 {MFD의 이름}/{PS의 이름}
으로 결정됩니다. (대소문자 구분)
이 글에서 예시로 든 것처럼 다음과 같이 SD를 설계하였다면, HC의 이름은 “Detector/EDep”이 되는 것입니다.
사용자는 이 HC의 이름만 있으면, G4SDManager로부터 GetCollectionID() 함수를 이용해 해당 HC의 고유 ID번호를 가져올 수 있습니다.
이 부분은 약간 C++ 측면의 깊은 설명이 필요합니다.
먼저 아무도 사용하지 않는 원초적인 방법부터 보여드리겠습니다.
1G4int fHCID = G4SDManager::GetSDMpointer()->GetCollectionID("Detector/EDep");
이렇게하면 fHCID
라는 변수에 "Detector/EDep"이라는 HC의 고유 ID번호가 저장됩니다. 하지만, Geant4에서는 아무도 이렇게 사용하지 않습니다.
GetCollectionID() 함수는 시간을 많이 잡아먹는 느린 함수입니다. 이런 함수를 매 Event가 종료될 때마다 실행하도록 하면 그만큼 시뮬레이션에 소요되는 시간이 길어지겠지요. 이런 문제를 해결하고자, 맨 처음에만 HC의 ID를 찾아오고 그 이후에는 전에 찾아온 값을 재활용하기 위해 다음과 같이 코드를 작성하여 이용합니다.
include/EventAction.hh 파일을 열고, G4int 자료형의 fHCID라는 이름의 멤버변수를 만듭니다.
멤버변수가 뭔지, 왜 hh파일이라는 곳에 쓰는지 잘 모르셔도 상관없습니다. 나중에 C++의 클래스 개념을 조금만 공부하시면 금방 이해될 것입니다.
일단 include 폴더에 있는 EventAction.hh 파일을 열어봅니다.
그리고 다음 그림과 같이 private:
라고 쓰여있는 줄과 };
가 적혀있는 줄 사이에, G4int fHCID;
라고 적어서 G4int 자료형의 fHCID 멤버변수를 선언합니다.
(이미 만들어져 있다면 그냥 두시면 됩니다. 강의자료 제작 과정에서 혼선이 있었던 관계로, 이미 적혀있는 코드를 받으셨을 수도 있습니다)
다시 src/EventAction.cc 파일로 돌아옵니다.
헤더를 포함시키는 부분(#include ...
) 바로 밑에 보면, EventAction 클래스의 생성자가 있을 것입니다. EventAction::EventAction() ...
과 같이 적혀있는 부분입니다.
여기에 다음 그림과 같이 , fHCID(-1)
이라는 문구를 추가해줍니다.
이는 fHCID라는 멤버변수를 -1로 초기화하겠다는 의미입니다.
EventAction.cc 파일의 맨 위에 G4SDManager.hh 헤더를 포함시킵니다.
이는 G4SDManager 클래스의 GetCollectionID 함수를 사용하기 위함입니다.
마지막으로 EndOfEventAction 함수에 다음 두 줄의 코드를 작성합니다.
이렇게 하면 fHCID의 값이 -1일 때만(맨 처음) GetCollectionID 함수를 이용하여 fHCID를 적절한 HC의 ID로 설정하게 되고, 이후에는 fHCID에 -1이 아닌 다른 값이 들어있으므로 더이상 if문 내부가 실행되지 않습니다.
1if(fHCID == -1)
2 fHCID = G4SDManager::GetSDMpointer()->GetCollectionID("Detector/EDep");
이어서, 이 ID를 이용하여 HCofThisEvent로부터 HC를 꺼내옵니다.
EventAction.cc 파일의 맨 위에 G4THitsMap.hh 헤더를 포함시켜 주시고, EndOfEventAction 함수 안쪽에 다음의 한 줄만 써주시면 됩니다.
1auto hitsMap = static_cast<G4THitsMap<G4double> *>(HCE->GetHC(fHCID));
이상한 문구가 너무 많이 있어서 생소하실 수 있겠지만, 해석해보면 간단합니다.
HCE->GetHC(fHCID)
HCE(HCofThisEvent)로부터 fHCID라는 ID를 가진 HC를 가져옴
static_cast<G4THitsMap
앞서 가져온 HC를 G4THitsMap<G4double> *
자료형으로 형변환함
auto hitsMap = static_cast<G4THitsMap
이걸 hitsMap이라는 변수명으로 저장
이제 모르는 내용은 G4THitsMap<G4double> *
이라는 처음보는 자료형 뿐이군요. 이는 HitsCollection의 종류 중 하나입니다. 이전 글에서 Hit을 저장하는 방법이 여러가지가 있지만, PS를 이용하는 경우에는 값을 누적하여 하나의 값으로 저장하며, 이를 Copy Number라는 사물함 번호로 분류하여 저장한다고 했었지요.
C++에서 이처럼 사물함 번호마다 값을 분류하여 저장하는 데이터 저장 방식의 대표적인 예로 Map이 있습니다. 그래서 Geant4에서는 이 Map이라는 구조를 활용하여 Hit을 저장하기 위한 Map이라는 이름으로 G4THitsMap이라는 클래스를 사용하고 있으며, 그 Map에 저장될 데이터가 G4double(실수)형이라는 것을 명시하기 위해 G4THitsMap<G4double>
과 같이 작성하게 됩니다. 끝의 *
는 이 자료형의 포인터형이라는 뜻입니다.
이 설명을 잘 모르겠으면 그냥 넘어가시고, 나중에 C++ 공부를 좀 더 한 뒤에 이해하셔도 됩니다. 핵심은 다음의 두 가지입니다. 이거만 기억하셔도 충분합니다.
1auto hitsMap = static_cast<G4THitsMap<G4double> *>(HCE->GetHC(fHCID));
라고 쓰면, hitsMap이라는 변수명으로 HitsCollection을 가져올 수 있음
이 hitsMap에는 Copy Number라는 사물함 번호마다, PS가 누적해서 기록한 값이 저장되어 있음
앞서 가져온 HitsCollection에는 PS가 기록한 값이 Copy Number라는 번호로 분류된 사물함마다 저장되어 있을 것입니다.
어떤 사용자는 Sensitive Detector를 딱 하나의 지오메트리에만 달았을 수도 있고, 어떤 사용자는 여러 개의 지오메트리에 달았을 수도 있습니다. 게다가, 이 Event에서 입자들이 그 지오메트리들을 아예 안지나갔을 수도 있고(Hit이 0개), 모두 다 지나갔을 수도 있습니다.
이런 이유때문에, PS를 이용하여 설계한 SD를 이용한 경우에는 HitsCollection 안에 몇 개의 Hit이 들어있을지 아무도 모릅니다. 이런 모든 경우를 다 아우르기 위해, HC안에 들어있는 Hit을 모두 다 훑으며 확인하는 방식을 이용합니다. 이는 C++에서 제공하는 for-each 반복문을 사용하면 쉽게 해결됩니다.
1for (const auto &iter : *(hitsMap->GetMap()))
2{
3 // ...
4}
위와 같이 작성하면, hitsMap 안에 있는 모든 사물함들을 iter라는 변수명으로 접근할 수 있게 됩니다. 이 때 iter.first
는 사물함 번호가 되고, *(iter.second)
는 그 사물함 안에 들어있는 hit의 누적 값이 됩니다.
예를 들어, G4PSEnergyDeposit을 통해 기록된 deposit된 에너지의 총합을 살펴보고, 이 값이 0보다 큰 경우 터미널 화면에 출력하고자 한다면 다음과 같이 코드를 작성합니다.
1for (const auto &iter : *(hitsMap->GetMap()))
2{
3 auto eDep = *(iter.second);
4 if (eDep > 0.)
5 {
6 G4cout << "--- Energy Deposit:" << eDep / MeV << G4endl;
7 }
8}
G4cout
과 G4endl
은 각각 std::cout
과 std::endl
과 동일하다고 생각하셔도 무방합니다.
아마 이렇게만 쓰시면 분명 MeV
라는게 무엇인지 모르겠다며 에러가 뜰 것입니다!
앞서 지오메트리 파트에서 단위를 사용하려면 G4SystemOfUnits.hh라는 헤더가 필요하다고 했었죠?
여기서도 MeV
를 사용해야 하므로, EventAction.cc 파일의 맨 위에 #include "G4SystemOfUnits.hh"
라는 내용을 추가해주셔야 합니다.
위의 코드를 해석하자면 다음과 같이 되겠군요
*(iter.second)
)을 eDep이라는 변수명으로 저장이런 방법을 통해 PS가 저장한 값을 Event가 끝날 때 가져와서 확인할 수 있게 됩니다. 지금까지 수행한 내용을 모두 작성하면 EndOfEventAction 함수 내부는 다음과 같은 모습일 것입니다.
이렇게 작성하였을 때 출력되는 결과를 살펴보겠습니다. 개수가 너무 적으면 확률상 원하는 결과가 나오지 않을수도 있으니, run.mac 파일에서 /run/beamOn 100
이라고 되어있는 줄을 /run/beamOn 1000
과 같이 약간 늘려주세요.
build 디렉토리에 들어와서, 다음 명령어를 통해 tracking verbose 1단계와 함께 출력해서 살펴봅시다.
1make
2./g4_minimal run.mac > vb.out
입자가 phantom에 들어가지 않은 다음과 같은 경우에는 추가적인 출력이 없는 것을 확인할 수 있습니다.
이제 phantom이라는 문구를 검색해 보겠습니다.
다음과 같이 phantom을 지나가긴 했지만, 별다른 반응이 일어나지 않아 deposit된 에너지가 없는 경우(eDep > 0.
조건을 만족하지 못함)에도 추가적인 출력이 없는 것을 확인할 수 있습니다.
좀 더 찾다보면 다음처럼 phantom에서 반응이 일어난 경우도 찾을 수 있습니다.
노란색으로 표시한 각 반응에서의 dE(MeV)들을 한 Event 내에서 다 합쳐져서, 그 총 합 값이 Event가 끝난 뒤에 --- Energy Deposit:...
와 같이 출력되고 있음을 확인할 수 있습니다.
이번 글에서 작성한 코드는 이 링크를 통해 다운받을 수 있습니다.
혹은 git repository를 clone하신 분의 경우에는, example branch의 이전 커밋 중 V3 scoring이라는 커밋을 참고하셔도 됩니다.
스코어링이 사실 쉽지 않습니다. 저도 많은 분들에게 강의를 수 차례 해보았지만, 제일 어려워하시는 부분이 스코어링입니다. C++의 문법을 모르면 생소한 코드가 너무 많아서 더욱 어려워보이는 것 같기도 합니다.
어떻게 하면 더 쉽게 전달할 수 있을까 많이 고민해 보았습니다만, 아직도 쉬워보이지는 않네요
하지만 스코어링 이론을 잘 숙지하고 있으시다면, 몇 번만 직접 구현해서 활용해보면 금방 터득하실 수 있으리라 생각합니다.
이제 다음 내용이 이 시리즈의 마지막입니다. 다음 글부터는 스코어링 결과를 파일로 출력하는 법에 대해 알아보겠습니다.