蒙特卡洛模拟GEANT4例子exampleB4a源码解读.doc

《蒙特卡洛模拟GEANT4例子exampleB4a源码解读.doc》由会员分享，可在线阅读，更多相关《蒙特卡洛模拟GEANT4例子exampleB4a源码解读.doc（29页珍藏版）》请在三一文库上搜索。

1、Geant4 ExampleB4a 源码解读今天继续简单讲一下exampleB4中的几个子例子，先是exampleB4a，仍然是采用Steppi ngAction 步一步抽取事件，这种方法的好处是能够很详细的考虑、输出每一步的信 息，但有时候我们只希望对感兴趣部分（位置）进行事件抽取，则用SD和Hit要明显方便一些，例如在exampleB4c中用到的就是 SD和Hit进行数据抽取。exampleB4a和exampleB4c 描述的是同一个几何、物理过程，只是数据抽取的方式不同，因此将这两个例子放在一起具 有比较意义。本次先对 exampleB4a进行讲解。exampleB4a计算了入射粒子（默

2、认电子，能 量50 MeV，轴向方向）在多层铅-液氩材料中（由10个吸收体和间隙复制组成）的能量沉 积及带电粒子（包括次级粒子）径迹长度。图1.几何可视化及粒子输运可视化exampleB4a.cc#in clude B4DetectorCo nstructio n.hh#in clude B4aActio nlni tializatio n.hh#ifdef G4MULTITHREADED#in clude G4MTRu nMan ager.hh#else#in clude G4Ru nMan ager.hh#en dif#in clude G4UIma nager.hh#in clude G

3、4UIcomma nd.hh#include FTFP_BERT.hh #include Randomize.hh#include G4VisExecutive.hh#include G4UIExecutive.hh/oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOooonamespace void PrintUsage() G4cerr Usage: G4endl;G4cerr exampleB4a -m macro -u UIsession -t nThreads G4endl;G4cerr note: -t option is available onl

4、y for multi-threaded mode. 7 ) PrintUsage(); return 1;G4String macro;G4String session;#ifdef G4MULTITHREADEDG4int nThreads = 0;#endiffor ( G4int i=1; i 0 ) runManager-SetNumberOfThreads(nThreads);#elseauto runManager = new G4RunManager;#endif/ Set mandatory initialization classes/auto detConstructio

5、n = new B4DetectorConstruction(); runManager-SetUserInitialization(detConstruction);auto physicsList = new FTFP_BERT; runManager-SetUserInitialization(physicsList);auto actionInitialization = new B4aActionInitialization(detConstruction); runManager-SetUserInitialization(actionInitialization);/ Initi

6、alize visualization/auto visManager = new G4VisExecutive;/ G4VisExecutive can take a verbosity argument - see /vis/verbose guidance./ G4VisManager* visManager = new G4VisExecutive(Quiet); visManager-Initialize();/ Get the pointer to the User Interface managerauto UImanager = G4UImanager:GetUIpointer

7、);/ Process macro or start UI session/if ( macro.size() ) / batch modeG4String command = /control/execute ;UImanager-ApplyCommand(command+macro);else / interactive mode : define UI sessionUImanager-ApplyCommand(/control/execute init_vis.mac);if (ui-IsGUI() UImanager-ApplyCommand(/control/execute gu

8、i.mac);ui-SessionStart();delete ui;/ Job termination/ Free the store: user actions, physics_list and detector_description are/ owned and deleted by the run manager, so they should not be deleted/ in the main() program !delete visManager;delete runManager;用户在终端输入字符串个数大于 7 则调用 PrintUsage 提示出错，并返回 1，跳出

9、 main()函数。通过for循环对输入字符串扫描识别，若程序名后字符串为-m”给macro赋值该字符串后一个字符串；若程序名后字符串为“ -u”给session赋值该字符串后一个字符串，否则调用PrintUsage提示出错，并返回 1跳出main()函数。定义一个 UI操作类对象，ui指针 初始化空指针，若用户没有调用 macro 控制文件，则 new 一个 ui 指针，分内存空间辟。初 始化探测器、物理过程以及用户行为，其中，用户行为与探测器构造有关。new 一个可视化管理器并对其初始化。定义一个 UI 管理类，通过判断条件 macro 是否为空，来决定执行交 互式界面模式还是批处理模式。

10、最后delete所有管理类，释放内存。B4DetectorConstruction.hh#ifndef B4DetectorConstruction_h #define B4DetectorConstruction_h 1#include G4VUserDetectorConstruction.hh#include globals.hhclass G4VPhysicalVolume;class G4GlobalMagFieldMessenger;class B4DetectorConstruction : public G4VUserDetectorConstructionpublic:B4D

11、etectorConstruction();virtual B4DetectorConstruction();public:virtual G4VPhysicalV olume* Construct();virtual void ConstructSDandField();/ get methods/const G4VPhysicalVolume* GetAbsorberPV() const; /const 修饰指针， 指针的内存空间 数据不能改变；函数重载const G4VPhysicalV olume* GetGapPV() const;private:/ methods/void Def

12、ineMaterials();G4VPhysicalV olume* DefineV olumes();/ data members/static G4ThreadLocal G4GlobalMagFieldMessenger* fMagFieldMessenger;/ magnetic field messengerG4VPhysicalV olume* fAbsorberPV; / the absorber physical volumeG4VPhysicalV olume* fGapPV; / the gap physical volumeG4bool fCheckOverlaps; /

13、 option to activate checking of volumes overlaps;/ inline functionsinline const G4VPhysicalV olume* B4DetectorConstruction:GetAbsorberPV() const return fAbsorberPV;inline const G4VPhysicalV olume* B4DetectorConstruction:GetGapPV() const return fGapPV;B4DetectorConstruction 继承于 G4VUserDetectorConstru

14、ction 基类， 声明构造函数、 析构 函数，Con struct。返回物理体。虚方法Con structSDa ndField()用于定义SD探测器以及空间内 的电磁场。 GetAbsorberPV() 和 GetGapPV() 为自定义的内联函数在类体外进行定义，分别返 回物理体 fAbsorberPV 和 fGapPV (私有数据成员) 。 G4GlobalMagFieldMessenger 类用于描 述电磁场的大小、方向等信息。B4DetectorConstruction.ccB4DetectorConstruction:B4DetectorConstruction(): G4VUs

15、erDetectorConstruction(),fAbsorberPV(nullptr),fGapPV(nullptr),fCheckOverlaps(true)/oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOoooB4DetectorConstruction:B4DetectorConstruction()/oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOoooG4VPhysicalVolume* B4DetectorConstruction:Construct() / Define materialsDefineMat

16、erials();/ Define volumesreturn DefineVolumes();/oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOooovoid B4DetectorConstruction:DefineMaterials()/ Lead material defined using NIST Managerauto nistManager = G4NistManager:Instance(); nistManager-FindOrBuildMaterial(G4_Pb);/ Liquid argon materialG4double a; /

17、 mass of a mole;G4double z; / z=mean number of protons;G4double density;new G4Material(liquidArgon, z=18., a= 39.95*g/mole, density= 1.390*g/cm3);/ The argon by NIST Manager is a gas with a different density/ Vacuumnew G4Material(Galactic, z=1., a=1.01*g/mole,density= universe_mean_density, kStateGa

18、s, 2.73*kelvin, 3.e-18*pascal);/ Print materialsG4cout *(G4Material:GetMaterialTable() G4endl;/oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOoooG4VPhysicalVolume* B4DetectorConstruction:DefineV olumes()/ Geometry parametersG4int nofLayers = 10;/ 定义有多少层 copynumberG4double absoThickness = 10.*mm;G4double g

19、apThickness = 5.*mm;G4double calorSizeXY = 10.*cm;/整体厚度auto layerThickness = absoThickness + gapThickness;/一层厚度auto calorThickness = nofLayers * layerThickness;auto worldSizeXY = 1.2 * calorSizeXY;auto worldSizeZ = 1.2 * calorThickness;/ Get materialsauto defaultMaterial = G4Material:GetMaterial(Gal

20、actic);auto absorberMaterial = G4Material:GetMaterial(G4_Pb);auto gapMaterial = G4Material:GetMaterial(liquidArgon);if ( ! defaultMaterial | ! absorberMaterial | ! gapMaterial ) G4ExceptionDescription msg;msg Cannot retrieve materials already defined.;G4Exception(B4DetectorConstruction:DefineV olume

21、s(), MyCode0001, FatalException, msg);/ World/auto worldS= new G4Box(World,/ its nameworldSizeXY/2, worldSizeXY/2, worldSizeZ/2); / its sizeauto worldLV= new G4LogicalV olume(worldS, defaultMaterial, World);auto worldPV= new G4PVPlacement(0,G4ThreeVector(), worldLV , World, 0, false, 0, fCheckOverla

22、ps);/ Calorimeter/auto calorimeterS/ its solid/ its material/ its name/ no rotation/ at (0,0,0)/ its logical volume/ its name/ its mother volume / no boolean operation/ copy number/ checking overlaps/ its name= new G4Box(Calorimeter,calorSizeXY/2, calorSizeXY/2, calorThickness/2); / its sizeauto cal

23、orLV= new G4LogicalV olume(calorimeterS, defaultMaterial, Calorimeter);new G4PVPlacement(0,G4ThreeVector(), calorLV, Calorimeter, worldLV , false, 0, fCheckOverlaps);/ its solid/ its material/ its name/ no rotation/ at (0,0,0)/ its logical volume/ its name/ its mother volume/ no boolean operation /

24、copy number/ checking overlaps/ Layer/auto layerS= new G4Box(Layer,/ its namecalorSizeXY/2, calorSizeXY/2, layerThickness/2); / its sizeauto layerLV= new G4LogicalV olume(layerS,/ its soliddefaultMaterial, / its material Layer);/ its namenew G4PVReplica(Layer,/ its namelayerLV,/ its logical volumeca

25、lorLV,/ its motherkZAxis,/ axis of replicationnofLayers,/ number of replicalayerThickness); / witdth of replica/ Absorber/auto absorberS= new G4Box(Abso,/ its namecalorSizeXY/2, calorSizeXY/2, absoThickness/2); / its sizeauto absorberLV= new G4LogicalV olume(absorberS, / its solid absorberMaterial,

26、/ its material Abso);/ its namefAbsorberPV= new G4PVPlacement(0,/ no rotationG4ThreeVector(0., 0., -gapThickness/2), / its position absorberLV,/ its logical volumeAbso, layerLV, false,0, fCheckOverlaps);/ its name/ its mother volume/ no boolean operation/ copy number/ checking overlaps / Gap/auto ga

27、pS= new G4Box(Gap,/ its namecalorSizeXY/2, calorSizeXY/2, gapThickness/2); / its sizeauto gapLV/ its solid/ its material/ its name= new G4LogicalV olume( gapS, gapMaterial, Gap);fGapPV= new G4PVPlacement(0,/ no rotationG4ThreeVector(0., 0., absoThickness/2), / its position gapLV,/ its logical volume

28、Gap,/ its namelayerLV,/ its mother volumefalse,/ no boolean operation0,/ copy numberfCheckOverlaps); / checking overlaps/ print parameters/G4cout G4endl G4endl The calorimeter is nofLayers layers of: absoThickness/mm mm of GetName() + gapThickness/mm mm of GetName() G4endl SetVisAttributes (G4VisAtt

29、ributes:GetInvisible();auto simpleBoxVisAtt= new G4VisAttributes(G4Colour(1.0,1.0,1.0); simpleBoxVisAtt-SetVisibility(true);calorLV-SetVisAttributes(simpleBoxVisAtt);/ Always return the physical World/return worldPV;/oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOooovoid B4DetectorConstruction:ConstructSD

30、andField()/ Create global magnetic field messenger./ Uniform magnetic field is then created automatically if/ the field value is not zero.G4ThreeVector fieldV alue;fMagFieldMessenger = new G4GlobalMagFieldMessenger(fieldValue); fMagFieldMessenger-SetVerboseLevel(1);/ Register the field messenger for

31、 deleting G4AutoDelete:Register(fMagFieldMessenger);在 con struct。函数中，定义了 Calorimeter ”、Layer ”、“ Abso ” 以及Gap ”四个逻辑 体，它们之间的关系是：构建一个大长方体Calorimeter ”，往里面复制并放置 10个规格相同的长方体称之为层Layer ”，再往Layer ”中分别放置 Abso ”以及Gap ”。在ConstructSDandField() 函数体里，我们可以定义 SD 探测器以及空间电磁场， 在本例 exampleB4a中我们只是定义了电磁场，可以通过run 1.mac中

32、命令行控制磁场属性：/globalField/setValue 0.2 0 0 tesla/设置均匀磁场，强度0.2特斯拉，方向沿 x轴exampleB4c 中，我们在 ConstructSDandField() 函数体里还定义了两个 SD 探测器，这样 做简化了后续数据抽取工作中对每一个 step 的位置判断，具体怎么操作到 B4c 例中再说。B4aActionInitialization.hh#ifndef B4aActionInitialization_h#define B4aActionInitialization_h 1#include G4VUserActionInitializa

33、tion.hhclass B4DetectorConstruction;/ Action initialization class./class B4aActionInitialization : public G4VUserActionInitializationpublic:B4aActionInitialization(B4DetectorConstruction*);virtual B4aActionInitialization();virtual void BuildForMaster() const;virtual void Build() const;private:B4Dete

34、ctorConstruction* fDetConstruction;#endif在构造函数中，用到了 B4DetectorConstruction 类的对象做形参，这是因为在后续判断 粒子当前位置过程中用到了探测器几何， 用于获取几何的物理体， 后面会说， 这里只需要知 道用户初始化过程用到了 B4DetectorConstruction 类。B4aActionInitialization.ccB4aActionInitialization:B4aActionInitialization(B4DetectorConstruction* detConstruction): G4VUserAct

35、ionInitialization(),fDetConstruction(detConstruction) /oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOoooB4aActionInitialization:B4aActionInitialization() /oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOooovoid B4aActionInitialization:BuildForMaster() constSetUserAction(new B4RunAction);/oooOO0OOooooooOO0OOoooo

36、ooOO0OOooooooOO0OOooovoid B4aActionInitialization:Build() const SetUserAction(new B4PrimaryGeneratorAction);SetUserAction(new B4RunAction);auto eventAction = new B4aEventAction;SetUserAction(eventAction);SetUserAction(new B4aSteppingAction(fDetConstruction,eventAction);/在 step 里面用到了 fDetConstruction

37、 变量将 detConstruction 赋值给 fDetConstruction ，另外，对粒子发射器、 RunAction 、EventAction 进行用户行为设置，其中， B4aSteppingAction 用到了探测器几何信息，一般是用于位置判 断。SetUserAction(new B4aSteppingAction(fDetConstruction,eventAction);B4PrimaryGeneratorAction.hhclass G4ParticleGun;class G4Event;/ The primary generator action class with p

38、article gum./ It defines a single particle which hits the calorimeter/ perpendicular to the input face. The type of the particle/ can be changed via the G4 build-in commands of G4ParticleGun class / (see the macros provided with this example).class B4PrimaryGeneratorAction : public G4VUserPrimaryGen

39、eratorAction public:B4PrimaryGeneratorAction();virtual B4PrimaryGeneratorAction();virtual void GeneratePrimaries(G4Event* event);/ set methodsvoid SetRandomFlag(G4bool value);private:G4ParticleGun* fParticleGun; / G4 particle gun;粒子发射器用的是 G4ParticleGun 类。B4PrimaryGeneratorAction.cc#include B4Primary

40、GeneratorAction.hh#include G4RunManager.hh#include G4LogicalV olumeStore.hh#include G4LogicalV olume.hh#include G4Box.hh#include G4Event.hh#include G4ParticleGun.hh#include G4ParticleTable.hh#include G4ParticleDefinition.hh#include G4SystemOfUnits.hh#include Randomize.hh/oooOO0OOooooooOO0OOooooooOO0

41、OOooooooOO0OOoooB4PrimaryGeneratorAction:B4PrimaryGeneratorAction(): G4VUserPrimaryGeneratorAction(),fParticleGun(nullptr)G4int nofParticles = 1;fParticleGun = new G4ParticleGun(nofParticles);/ default particle kinematic/auto particleDefinition= G4ParticleTable:GetParticleTable()-FindParticle(e-); fParticleGun-SetParticleDefinition(particleDefinition); fParticleGun-SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.); fParticleGun-SetParticleEnergy(50.*MeV);/oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOoooB4PrimaryGeneratorAction:B4PrimaryGeneratorAction()