Geant4光学过程教程

基本原理

参考例子:extended/optices/OpNovice,建议新的探测器在此程序的基础上修改几何、材质和传感器。修改几何只需要修改GDML文件,基本不需要改动OpNovice;采样记录需要在恰当的过程中增加记录或者输出代码。

光学光子opticalphoton和gamma不同,在G4OpticalPhoton中描述。

光学过程包含以下头文件:

  • G4Cerenkov 产生切伦科夫光
  • G4Scintillation 产生闪烁光
  • G4OpAbsorption 吸收
  • G4OpRayleigh 散射
  • G4OpBoundaryProcess 界面过程,如反射折射

GDML基本结构

GDML的基本用法见 https://gdml.web.cern.ch/GDML/,参考`extended/optics/OpNovice`的自带gdml和`persistency/gdml/G01/optical

  • <define> </define>定义常数和变量
  • <materials> </materials>定义材料
  • <solids> </solids>定义solid
  • <structure> </structure>使用solid建立volume,并指定内部包含volume的位置
  • <setup> </setup>:指定World对应的volume

定义material的光学性质

材料的光学性质数据在<define/>中定义,见光学过程文档中的表格。比较常用的光学性质主要是RINDEX,材料的折射率,是一系列光子能量-折射率数值的数组。定义折射率后界面处的反射率会根据菲涅耳定律计算,也会考虑光学光子的偏振。界面有特殊需要的设定参见光学过程文档。

<material/>中附加在材料上。例如水:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
<material name="Water" state="solid">
<property name="RINDEX" ref="WATERRINDEX"/>
<property name="MIEHG" ref="WATERMIEHG"/>
<property name="ABSLENGTH" ref="WATERABSLENGTH"/>
<property name="SCINTILLATIONCOMPONENT1" ref="WATERSCINTILLATIONCOMPONENT1"/>
<property name="SCINTILLATIONCOMPONENT2" ref="WATERSCINTILLATIONCOMPONENT2"/>
<property name="MIEHG_FORWARD" ref="WATERMIEHG_FORWARD"/>
<property name="MIEHG_BACKWARD" ref="WATERMIEHG_BACKWARD"/>
<property name="MIEHG_FORWARD_RATIO" ref="WATERMIEHG_FORWARD_RATIO"/>
<property name="SCINTILLATIONYIELD" ref="WATERSCINTILLATIONYIELD"/>
<property name="RESOLUTIONSCALE" ref="WATERRESOLUTIONSCALE"/>
<property name="SCINTILLATIONTIMECONSTANT1" ref="WATERSCINTILLATIONTIMECONSTANT1"/>
<property name="SCINTILLATIONTIMECONSTANT2" ref="WATERSCINTILLATIONTIMECONSTANT2"/>
<property name="SCINTILLATIONYIELD1" ref="SCINTILLATIONYIELD1"/>
<property name="SCINTILLATIONYIELD2" ref="SCINTILLATIONYIELD2"/>
<T unit="K" value="293.15"/>
<MEE unit="eV" value="68.9984174679527"/>
<D unit="g/cm3" value="1"/>
<fraction n="0.112097669256382" ref="Hydrogen"/>
<fraction n="0.887902330743618" ref="Oxygen"/>
</material>

其中MEE指的是平均激发能量,和光学过程关系不大。温度可以不设定,默认为273.15K,密度必须指定。

Geant4光学过程在材料中产生切伦科夫光的光子能量范围是材料RINDEX的范围,在此范围外的切伦科夫光子不会产生。

定义surface的光学性质opticalsurface

表面光学性质在<solid/>中定义:

1
2
3
4
<opticalsurface name="surf1" model="glisur" finish="polished" type="dielectric_dielectric" value="1.0"/>
<opticalsurface name="surf2" model="glisur" finish="polished" type="dielectric_metal" value="1.0">
<property name="REFLECTIVITY" ref="REFLECTIVITY" />
</opticalsurface>

model是表面模型方案,默认为glisur;finish为表面平整度,默认为polished镜面反射,可以选漫反射ground等;type为表面两边的介电类型默认为dielectric-dielectric。可以选择的范围见源文件/persistency/gdml/src/G4GDMLReadSolids.cc。默认的设置在源文件/processes/optical/src/G4OpBoundaryProcess.cc中。

<structure/>中将opticalsurface附在volume上,其中skinsurface表示整个volume被该表面覆盖,bordersurface表示一对有序的volume之间的表面由这种表面决定。

1
2
3
4
5
6
7
<skinsurface name="AirSurface" surfaceproperty="AirSurface">
<volumeref ref="Bubble_log"/>
</skinsurface>
<bordersurface name="WaterSurface" surfaceproperty="WaterSurface">
<physvolref ref="Tank_phys"/>
<physvolref ref="expHall_phys"/>
</bordersurface>

为了证明对材料性质的定义确实产生了正确的光学过程,我们考察以下示例。

示例 单一截面折反射

在Geant4官方案例extended/optical/OpNovice基础上设计。利用这个案例我们设计一个最简单的一个光子穿过一块材料发生折射和反射的情景。

几何

探测器几何为空气中放置一个长方体介质,材料可以变化,先设成水。探测器GDML文件SimpleReflector.gdml如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
<?xml version="1.0" encoding="UTF-8" ?>
<gdml xmlns:gdml="http://cern.ch/2001/Schemas/GDML" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:noNamespaceSchemaLocation="http://service-spi.web.cern.ch/service-spi/app/releases/GDML/schema/gdml.xsd" >

<define>
<matrix coldim="2" name="AIRRINDEX" values="2.034*eV 1 2.068*eV 1 2.103*eV 1 2.139*eV 1 2.177*eV 1 2.216*eV 1 2.256*eV 1 2.298*eV 1 2.341*eV 1 2.386*eV 1 2.433*eV 1 2.481*eV 1 2.532*eV 1 2.585*eV 1 2.64*eV 1 2.697*eV 1 2.757*eV 1 2.82*eV 1 2.885*eV 1 2.954*eV 1 3.026*eV 1 3.102*eV 1 3.181*eV 1 3.265*eV 1 3.353*eV 1 3.446*eV 1 3.545*eV 1 3.649*eV 1 3.76*eV 1 3.877*eV 1 4.002*eV 1 4.136*eV 1"/>
<matrix coldim="2" name="WATERRINDEX" values="2.034*eV 1.3435 2.068*eV 1.344 2.103*eV 1.3445 2.139*eV 1.345 2.177*eV 1.3455 2.216*eV 1.346 2.256*eV 1.3465 2.298*eV 1.347 2.341*eV 1.3475 2.386*eV 1.348 2.433*eV 1.3485 2.481*eV 1.3492 2.532*eV 1.35 2.585*eV 1.3505 2.64*eV 1.351 2.697*eV 1.3518 2.757*eV 1.3522 2.82*eV 1.353 2.885*eV 1.3535 2.954*eV 1.354 3.026*eV 1.3545 3.102*eV 1.355 3.181*eV 1.3555 3.265*eV 1.356 3.353*eV 1.3568 3.446*eV 1.3572 3.545*eV 1.358 3.649*eV 1.3585 3.76*eV 1.359 3.877*eV 1.3595 4.002*eV 1.36 4.136*eV 1.3608"/>
</define>

<materials>
<isotope N="14" Z="7" name="N14">
<atom unit="g/mole" value="14.0031"/>
</isotope>
<isotope N="15" Z="7" name="N15">
<atom unit="g/mole" value="15.0001"/>
</isotope>
<element name="Nitrogen">
<fraction n="0.99632" ref="N14"/>
<fraction n="0.00368" ref="N15"/>
</element>
<isotope N="16" Z="8" name="O16">
<atom unit="g/mole" value="15.9949"/>
</isotope>
<isotope N="17" Z="8" name="O17">
<atom unit="g/mole" value="16.9991"/>
</isotope>
<isotope N="18" Z="8" name="O18">
<atom unit="g/mole" value="17.9992"/>
</isotope>
<element name="Oxygen">
<fraction n="0.99757" ref="O16"/>
<fraction n="0.00038" ref="O17"/>
<fraction n="0.00205" ref="O18"/>
</element>
<material name="Air" state="gas">
<property name="RINDEX" ref="AIRRINDEX"/>
<D unit="g/cm3" value="0.00129"/>
<fraction n="0.7" ref="Nitrogen"/>
<fraction n="0.3" ref="Oxygen"/>
</material>
<isotope N="1" Z="1" name="H1">
<atom unit="g/mole" value="1.00782503081372"/>
</isotope>
<isotope N="2" Z="1" name="H2">
<atom unit="g/mole" value="2.01410199966617"/>
</isotope>
<element name="Hydrogen">
<fraction n="0.999885" ref="H1"/>
<fraction n="0.000115" ref="H2"/>
</element>
<material name="Water" state="solid">
<property name="RINDEX" ref="WATERRINDEX"/>
<D unit="g/cm3" value="1"/>
<fraction n="0.112097669256382" ref="Hydrogen"/>
<fraction n="0.887902330743618" ref="Oxygen"/>
</material>

</materials>

<solids>
<box lunit="mm" name="World" x="100" y="100" z="100"/>
<box lunit="mm" name="Bulk" x="20" y="100" z="100"/>
</solids>

<structure>
<volume name="Bulk">
<materialref ref="Water"/>
<solidref ref="Bulk"/>
</volume>
<volume name="World">
<materialref ref="Air"/>
<solidref ref="World"/>
<physvol name="Bulk">
<volumeref ref="Bulk"/>
<position name="Bulk_pos" unit="mm" x="0" y="0" z="0"/>
</physvol>
</volume>
</structure>

<setup name="Default" version="1.0">
<world ref="World"/>
</setup>

</gdml>

这个GDML文件只定义了必要的光学参数,包括两种介质的折射率。由于没有规定界面性质,默认为glisur模型镜面、电介质—电介质。

为了便于验证折射定律,在src/OpNoviceSteppingAction.ccvoid OpNoviceSteppingAction::UserSteppingAction(const G4Step* step)if块的最后一行增加输出每次粒子穿过界面进入新的volume时的动量方向矢量:

1
G4cout<<"Next Direction:"<<step->GetTrack()->GetDynamicParticle()->GetMomentumDirection().x()<<" "<<step->GetTrack()->GetDynamicParticle()->GetMomentumDirection().y()<< " "<<step->GetTrack()->GetDynamicParticle()->GetMomentumDirection().z()<<G4endl;

编译该例子。

为了方便运行测试,写一个宏文件op.mac,设定入射粒子为光学光子,能量为2.034eV(609.5nm,定义的折射率为1.3435),方向在xy平面内。

1
2
3
4
/gun/particle opticalphoton
/gun/energy 2.034 eV
/gun/position -40 0 0 mm
/gun/direction 1 0.5 0

运行例子:

1
./OpNovice -g ../gdml/SimpleReflector.gdml

在Session中载入宏文件:

1
/control/execute op.mac

可以通过打开轨迹输出来观察粒子运行轨迹(可以增加输出的小数点位数):

1
/tracking/verbose 1

可以通过打开光学边界过程输出来观察粒子和边界的交互:(0 = silent; 1 = initialisation; 2 = during tracking)

1
/process/optical/boundary/verbose 1

点击运行按钮或者在Session中输入:

1
/run/beamOn 1

沿指定方向产生一个光子。

在输出中,得到第一个折射方向数据为:

1
> Next Direction:0.942972003905678 0.3328720472645761 0

即在xy平面内,入射动量方向(1,0.5),折射动量方向(0.942972003905678,0.3328720472645761),可以通过斯涅耳公式求出折射率为0.332872,恰为我们在GDML文件中设定的值。改变光学光子能量,能够算出对应的折射率。

在这样的设定下,Geant4除了产生折射的轨迹,也有概率产生反射的轨迹,考虑到只要折射率确定了,介质界面处的透射率就确定了(菲涅耳定律)。对此,可以在Session中输入:

1
/run/beamOn 10

产生多个光学光子,产生的轨迹为