Geant4 Tips: Optical/Scintillation

HardRyan

已于 2024-04-06 20:47:14 修改

阅读量881

点赞数 12

分类专栏： Geant4 文章标签：笔记

于 2024-04-06 20:44:22 首次发布

本文链接：https://blog.csdn.net/m0_51809467/article/details/137402124

版权

Geant4 专栏收录该内容

12 篇文章 3 订阅

订阅专栏

写在前面

出于需要，在此整理学习Geant4 G4Scintillation的笔记，以及整理部分G4Optical Physical的相关内容。笔者作为工科材料专业，加上年级较低能力不足，部分理论无法完全解释，因此对于本次笔记的内容安排如下：

以个人感兴趣的内容为主，即G4Scintillation, G4OpticalSurface
适当罗列其他同在 G4OpticalPhysicsList相关的内容，仅整理最基本的用法

笔记内容来源于官网以及\examples\extended\optical\。

G4OpticalPhysicsList

在介绍G4Scintillation前，先简单介绍以下G4OpticalPhysicsList的使用，首先要注册，如下：

#include "G4OpticalPhysics.hh"
G4VModularPhysicsList* physicsList = new FTFP_BERT;  // for example
G4OpticalPhysics* opticalPhysics = new G4OpticalPhysics();
physicsList->RegisterPhysics(opticalPhysics);
runManager->SetUserInitialization(physicsList);

然后需要声明一些Optical Process 需要的参数，相应的参数如果存在，则可能调用这些过程，反之则不会。

例如为了使用Cerenkov effect，需要定义材料的折射率，则以下代码在可能的情况下可以激活Cerenkov effect，

G4MaterialPropertiesTable* myMPT1 = new G4MaterialPropertiesTable();
myMPT1->AddProperty("RINDEX", &photonEnergy[0], &refractiveIndex1[0],numEntries);
water->SetMaterialPropertiesTable(myMPT1);

AddProperty

在上一段代码中，可能唯一需要花时间理解的代码就是第二行了，事实上 AddProperty 还有一个很像的”兄弟“ AddConstantProperty。

在 include 中，这两个方法的原型是：


  void AddConstProperty(const G4String& key, G4double propertyValue,
                        G4bool createNewKey = false);

  G4MaterialPropertyVector* AddProperty(
    const G4String& key, const std::vector<G4double>& photonEnergies,
    const std::vector<G4double>& propertyValues, G4bool createNewKey = false,
    G4bool spline = false);
  // Add a new property to the table by giving a key-name and
  // vectors of values

显然还有其他方式，例如在上一节中使用的是

  G4MaterialPropertyVector* AddProperty(
    const char* key, G4double* photonEnergies, G4double* propertyValues,
    G4int numEntries, G4bool createNewKey = false, G4bool spline = false);
  // Add a new property to the table by giving a key-name and the
  // arrays x and y of size NumEntries.

这种方式需要通过AddEntry进一步填充其中的键值对。

AddConstantProperty 和 AddProperty 从代码的注释中也可以看出，两者的主要区别是，前者的值不依赖于能量，或者说是光的波长，而后者依赖，所以需要输入两个对应的 std::<G4double> 。这种方式也比先AddProperty再AddEntry 来的方便和符合直觉。需要注意的是，两种方式还有一处不同在于后者（AddEntry）产生的属性默认spline（插值）和 createNewKey 都是 false。

使用的大概模板是：

std::vector<G4double> photonEnergy = {0., 1., 2., 3.};
std::vector<G4double> absorption = {0., 1., 2., 3. };
G4MaterialPropertiesTable* myMPT1 = new G4MaterialPropertiesTable();
myMPT1->AddProperty("ABSLENGTH", photonEnergy, absorption, false, true); 
myMPT1->AddProperty("SCINTILLATIONYIELD", 50./MeV);

// Setting Other Propertyes
// ...

myMaterial->SetMaterialPropertiesTable(myMPT1);

G4Scintillation

在进一步介绍其中各个 Process 需要的材料参数的设置前，需要说明前文出现的 “RINDEX”, “ABSLENGTH”, “SCINTILLATIONYIELD”，这些是与Process相关的参数设置的关键词，使用者也可以自定义关键词，并与自己创造的一些 Process 相互关联。

在Geant4中，闪烁体可以用单位能量的光子产额，产生的光子的能量分布，以及光子衰减的时间系数来描述。在此基础上还有两种特殊的应用情况，首先是多阶段的衰减，其次是粒子选择性的闪烁行为。这些将在后文分别说明。

SCINTILLATIONYIELD

闪烁体的光子产额与由关键词 "SCINTILLATIONYIELD” 定义的材料参数有关。Step产生的平均光子数是这个参数与其他两个参数的沉积，这两个参数分别是比例系数和能量，能量沉积的总能量，如果需要拟合非线性响应，则需要声明，示例如下

fLXe->GetIonisation()->SetBirksConstant(0.126*mm/MeV);

则使用的能量为VisibleEnergyDepositAtAStep 。

代码中的常数为Birks常数，对应与下列代码中的 $B$
$\frac{\mathrm{d}L}{\mathrm{d} x} = S \frac{\frac{\mathrm{d} E}{\mathrm{d}x}}{1+kB\frac{\mathrm{d} E}{\mathrm{d}x}}$
$L$ 为光产额， $S$ 为闪烁效率， $\frac{\mathrm{d} E}{\mathrm{d}x}$ 为单位长度辐射粒子乘积s的能量。

为了使最终平均光子数为预设的值，如果实际产生的光子数小于10，将会以通过由同样均值的泊松分布决定最终使用的光子数，如果大于10，则会通过相同均值、方差为均值的根与分辨率的乘积的高斯分布决定。r分辨率由关键词 “RESOLUTIONSCALE” 设置。

SCINTILLATIONYIELD 设置的示例为

myMPT1->AddConstantProperty("SCINTILLATIONYIELD", 50./MeV);
myMPT1->AddConstantProperty("RESOLUTIONSCALE", 0;

SCINTILLATIONCOMPONENT

用于设置闪烁光的光谱分布，需要获得光子的能谱分布。按照不同的衰减速度可以将发射光分为不同的组成部分，不同的组成部分可能能谱分布、衰减速度、光产额都不同

设置发射光的能谱的示例：

  G4double scintilFastArray[]{ 1.0, 1.0 };
  G4double energyArray[]{ 2.034 * eV, 4.136 * eV };
  G4int lenArray = 2;

  std::vector<G4double> scintilSlow = {
    0.01, 1.00, 2.00, 3.00, 4.00, 5.00, 6.00, 7.00, 8.00, 9.00, 8.00,
    7.00, 6.00, 4.00, 3.00, 2.00, 1.00, 0.01, 1.00, 2.00, 3.00, 4.00,
    5.00, 6.00, 7.00, 8.00, 9.00, 8.00, 7.00, 6.00, 5.00, 4.00
  };

  std::vector<G4double> photonEnergy = {
    2.034 * eV, 2.068 * eV, 2.103 * eV, 2.139 * eV, 2.177 * eV, 2.216 * eV,
    2.256 * eV, 2.298 * eV, 2.341 * eV, 2.386 * eV, 2.433 * eV, 2.481 * eV,
    2.532 * eV, 2.585 * eV, 2.640 * eV, 2.697 * eV, 2.757 * eV, 2.820 * eV,
    2.885 * eV, 2.954 * eV, 3.026 * eV, 3.102 * eV, 3.181 * eV, 3.265 * eV,
    3.353 * eV, 3.446 * eV, 3.545 * eV, 3.649 * eV, 3.760 * eV, 3.877 * eV,
    4.002 * eV, 4.136 * eV
  };
  myMPT1->AddProperty("SCINTILLATIONCOMPONENT1", energyArray, scintilFastArray,lenArray);
  myMPT1->AddProperty("SCINTILLATIONCOMPONENT2", photonEnergy, scintilSlow,false, true);
  myMPT1->AddConstProperty("SCINTILLATIONTIMECONSTANT1", 1. * ns);
  myMPT1->AddConstProperty("SCINTILLATIONTIMECONSTANT2", 10. * ns);
  myMPT1->AddConstProperty("SCINTILLATIONYIELD1", 0.8);
  myMPT1->AddConstProperty("SCINTILLATIONYIELD2", 0.2);

上面的示例不仅展示了不同的成分的发射光能谱的设置，还展示了相应的衰减时间和光产额的占比的设置。

Particle type

在某些应用情景，可能存在多种激发粒子，不同的激发粒子产生的闪烁行为也不同。设置相当简单，在相应的参数前加上粒子名称即可，例如

  myMPT1->AddConstProperty("ALPHASCINTILLATIONYIELD" , 10./MeV);

如果设置了相应的粒子对应的闪烁行为，则会优先使用粒子对应的参数。否则会使用通用的参数，例如，GAMMASCINTILLATIONYIELD没有设置，则会使用 SCINTILLATIONYIELD作为gamma粒子的光产额。

其他

闪烁过程还有一些其他可以设定的参数

Enable particle-dependent yields
Record track information
Whether to track secondaries before resuming tracking of primary particle
Whether to use a finite rise time for the secondary emission time
Whether to add produced optical photons to the stack (the alternative is to kill them)
Set the verbosity level during tracking

G4OpticalSurface

有“三”种能够调用特定物理过程的界面的方法。分别是通过设置"RINDEX"自动声明，和 G4LogicalSkinSurface 和 G4LogicalBorderSurface。

对于第一种方法，只需要设置所有相关volume的折射系数，Geant4会自动考虑在界面出的折射发射等现象。

对于 G4LogicalSkinSurface ，是声明包围一个volume的闭合界面为Opticalsurface，比较适应于表面改性或者涂覆光学材料等的场景，示例代码：

G4OpticalSurface* opAirSurface = new G4OpticalSurface("AirSurface");
auto bubbleAir_log = new G4LogicalVolume(bubbleAir_box, air, "Bubble", 0, 0, 0);
G4LogicalSkinSurface* airSurface = new G4LogicalSkinSurface("AirSurface", bubbleAir_log, opAirSurface);

对于 G4LogicalBorderSurface, 声明两个Physical Volume 接触的界面为光学界面，因为两个Physcial Volune 在声明的时候存在顺序，因此可以设置从Volume A到Volume B 和 B到A触发不同的光学过程，示例代码

auto expHall_phys = new G4PVPlacement(0, G4ThreeVector(), expHall_log, "World", 0, false, 0);
auto waterTank_phys = new G4PVPlacement( 0, G4ThreeVector(), waterTank_log, "Tank", expHall_log, false, 0);
G4OpticalSurface* opWaterSurface = new G4OpticalSurface("WaterSurface");
G4LogicalBorderSurface* waterSurface = new G4LogicalBorderSurface( "WaterSurface", waterTank_phys, expHall_phys, opWaterSurface);

最简单的界面过程的模拟只需要设置好“RINDEX”即可，如果需要更深的处理，才需要设置对应的OpticalSurface，选择对应的模型，设置对应的MaterialPropertiesTable等。

Model, Finish and type

在声明一个 G4OpticalSurface 后，可以选择不同的模型描述表面的类型，例如不同材料不同粗糙程度的计算模型，对应于 SetModel，在Model中，可以选择不同的光滑程度，对应于 SetFinish。设置界面的材料类型, 对应于 SetType。

SetModel 和 SetModel 的目的都是通过引入额外的信息数据，来尽量地使得拟合结果接近现实数据。如果只是做一些不深入的模拟过程，SetType很多时候就能解决问题。

以下是示例代码：

auto opAirSurface = new G4OpticalSurface("AirSurface");
opAirSurface->SetType(dielectric_dielectric);
opAirSurface->SetFinish(polished);
opAirSurface->SetModel(glisur);

模型的选择，和模型下界面以及其他参数的选择可见官网。

如果不设置模型，程序对界面的处理就是绝对光滑、两面均为绝缘体的界面。

其他过程

以下过程做简单罗列

Cerenkov Effect

当粒子在介质中的传播速度大于光在该介质传播的群速度时，就会产生 Cerenkov 辐射。

要让程序考虑 Cerenkov 过程，必须要要通过"RINDEX"标注材料的折射系数。

Cerenkov 过程还有许多其他能调整的参数。见官网。

Set the step size to limit the number of photons produced (on average) to a given value
Set the maximum change in β=v/c in a step
Specify whether to add Cerenkov photons to the stack, and track them.
Specify whether to track secondaries produced in the step before continuing with primary.
Set the verbosity of the process during tracking

Absorption

光子的吸收由平均自由程定义，是一个光子在被吸收前运动路径长度的平均距离，被关键词“”ABSLENGTH" 定义。

Rayleigh Scattering

瑞利散射过程需要定义相应的衰减距离，定义为光子在被瑞利散射前运动的平均距离，关键词为"RAYLEIGH"。如果没有提供瑞利散射对应的距离，部分情况下也可以用Geant4自动计算得出，此时需要设置"ISOTHERMAL_COMPRESSIBILITY" 参数或者设置材料为"Water"。最终的自由程还可以通过“RS_SCALE_FACTOR”比例系数进行调整。

Wavelength Shifting

WLS光纤是一种特殊的光纤，能够吸收某一波长的光线并发射为另一波长的光线。将WLS光纤耦合在光路在光路中可以抑制探测器的重吸收现象，提高光电倍增管受到的信号。

WLS过程可以由三个参数描述，分别式吸收、发射光谱和其中的时间延迟，对应于Geant4中的关键词分别是"WLSABSLENGTH"，“WLSCOMPONENT", “WLSTIMECONSTANT”

Set the time profile of the first WLS process to be either “delta” or “exponential”:
Set the verbosity of the first WLS process. 0 = silent; 1 = initialisation; 2 = during tracking.
Set the time profile of the second WLS process to be either “delta” or “exponential”.
Set the verbosity of the second WLS process. 0 = silent; 1 = initialisation; 2 = during tracking.

Mie Scattering

Mie Scattering 是用于球形颗粒对光子散射的麦克斯韦公式的分析解，使用范围是球形颗粒的尺度与光子的波长在相同的量级。需要提供相应的参数"MIEHG", “MIEHG_FORWARD”, “MIEHG_BACKWARD”, “MIEHG_FORWARD_RATIO”

总结

一个对以上提及对所有过程的代码示例：

G4MaterialPropertiesTable* myMPT1 = new G4MaterialPropertiesTable();

// refractive index
std::vector<G4double> photonEnergy = {2.034 * eV, 2.068 * eV, 2.103 * eV};
std::vector<G4double> refractiveIndex = {1.3435, 1.344,  1.3445};
myMPT1->AddProperty("RINDEX", photonEnergy, refractiveIndex);

// absorption
std::vector<G4double> absorption = {3.448 * m,  4.082 * m,  6.329 * m};
myMPT1->AddProperty("ABSLENGTH", photonEnergy, absorption);

// scintillation
// scintillation component
std::vector<G4double> scintilFastArray = { 1.0, 1.0 };
std::vector<G4double> energyArray = { 2.034 * eV, 4.136 * eV };
std::vector<G4double> scintilSlow = {0.01, 1.00, 2.00, 3.00};
myMPT1->AddProperty("SCINTILLATIONCOMPONENT1", energyArray, scintilFastArray,false, false);  // createNewKey: false, spline: false
myMPT1->AddProperty("SCINTILLATIONCOMPONENT2", photonEnergy, scintilSlow, false, true);

// scintillation constants
myMPT1->AddConstProperty("SCINTILLATIONYIELD", 50. / MeV);
myMPT1->AddConstProperty("RESOLUTIONSCALE", 1.0);
myMPT1->AddConstProperty("SCINTILLATIONTIMECONSTANT1", 1. * ns);
myMPT1->AddConstProperty("SCINTILLATIONTIMECONSTANT2", 10. * ns);
myMPT1->AddConstProperty("SCINTILLATIONYIELD1", 0.8);
myMPT1->AddConstProperty("SCINTILLATIONYIELD2", 0.2);

myMaterial->SetMaterialPropertiesTable(myMPT1);

// Mie
// Mie energy-dependent
std::vector<G4double> energy_water = {1.56962 * eV, 1.58974 * eV, 1.61039 * eV};
std::vector<G4double> mie_water = {167024.4 * m, 158726.7 * m, 150742 * m};
myMPT1->AddProperty("MIEHG", energy_water, mie_water, false, true);

// Mie constants
G4double mie_water_const[3] = { 0.99, 0.99, 0.8 };
myMPT1->AddConstProperty("MIEHG_FORWARD", mie_water_const[0]);
myMPT1->AddConstProperty("MIEHG_BACKWARD", mie_water_const[1]);
myMPT1->AddConstProperty("MIEHG_FORWARD_RATIO", mie_water_const[2]);


myMaterial->SetMaterialPropertiesTable(myMPT1);

// Set the Birks Constant for the scintillator
myMaterial->GetIonisation()->SetBirksConstant(0.126 * mm / MeV);

// optical surface
// optical skin surface
G4OpticalSurface* opSurface = new G4OpticalSurface("opSurface");
G4LogicalSkinSurface* Surface = new G4LogicalSkinSurface("skinSurface", theLogicalVolume, opSurface);

// optical border surfae
G4OpticalSurface* opSurface2 = new G4OpticalSurface("opSurface2");
G4LogicalBorderSurface* Surface2 = new G4LogicalBorderSurface("borderSurface", thePhysicalVolume1, thePhysicalVolume2, opSurface2);