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

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

2、子（默认电子，能 量50 MeV ,轴向方向）在多层铅-液僦材料中（由10个吸收体和间隙复制组成）的能量沉 积及带电粒子（包括次级粒子）径迹长度。图1.几何可视化及粒子输运可视化exampleB4a.ec#include B4DetectorConstruction.hh#include B4aActionInitialization.hh#ifdef G4MULTITHREADED#include G4MTRunManager.hh#else#include G4RunManager.hh#endif#include G4UImanager.hh#include G4UIcommand.hh

3、#include FTFP_BERT.hh”#include Randomize.hh#include G4VisExecutive.hh#include G4UIExecutive.hh.oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOspace void PrintUsage() G4cerr Usage: G4endl;G4cerr exampleB4a -m macro -u UIsession -t nThreads G4endl;G4cerr note: -t option is available only for multi-t

4、hreaded 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 detConstruction = new B4Det

5、ectorConstruction();runManager-SetUserInitialization(detConstruction);auto physicsList = new FTFP_BERT;runManager-SetUserInitialization(physicsList);auto actioninitialization = new B4aActionInitialization(detConstruction);runManager-SetUserInitialization(actionInitialization);/ Initialize visualizat

6、ion/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();/ Process macr

7、o 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 gui.mac);)ui-Sessi

8、onStart();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,跳出main()-m”给macro赋值该字

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

10、#ifndef B4DetectorConstruction_h#define B4DetectorConstruction_h 1#include G4VUserDetectorConstruction.hh#include globals.hhclass G4VPhysicalV olume;class G4GlobalMagFieldMessenger;class B4DetectorConstruction : public G4VUserDetectorConstruction public:B4DetectorConstruction();virtual B4DetectorCon

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

12、 olumes();/ data members/static G4ThreadLocal G4GlobalMagFieldMessenger* fMagFieldMessenger;/ magnetic field messengerG4VPhysicalV olume* fAbsorberPV; / the absorber physical volumeG4VPhysicalV olume* fGapPV; / the gap physical volumeG4bool fCheckOverlaps; / option to activate checking of volumes ov

13、erlaps ;/ inline functionsinline const G4VPhysicalV olume* B4DetectorConstruction:GetAbsorberPV() const return fAbsorberPV;inline const G4VPhysicalV olume* B4DetectorConstruction:GetGapPV() const return fGapPV;B4DetectorConstruction 继承于 G4VUserDetectorConstruction 基类，声明构造函数、析构 函数，Construct。返回物理体。虚方法

14、ConstructSDandField()用于定义SD探测器以及空间内 的电磁场。GetAbsorberPV()和GetGapPV()为自定义的内联函数在类体外进行定义，分别返 回物理体fAbsorberPV和fGapPV (私有数据成员)。G4GlobalMagFieldMessenger类用于描 述电磁场的大小、方向等信息。B4DetectorConstruction.ccB4DetectorConstruction:B4DetectorConstruction():G4VUserDetectorConstruction(),fAbsorberPV(nullptr),fGapPV(null

15、ptr), fCheckOverlaps(true)/.oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOooo.B4DetectorConstruction二B4DetectorConstruction()()/.oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOooo.G4VPhysicalVolume* B4DetectorConstruction:Construct()(/ Define materials DefineMaterials();/ Define volumes return DefineVolumes();

16、)/oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOooovoid B4DetectorConstruction:DefineMaterials()(/ Lead material defined using NIST Managerauto nistManager = G4NistManager二Instance。； nistManager-FindOrBuildMaterial(G4_Pb);/ Liquid argon materialG4double a; / mass of a mole;G4double z; / z=mean number of

17、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,kStateGas, 2.73*kelvin, 3.e-18*pascal);/ Print materials

18、G4cout *(G4Material:GetMaterialTable() G4endl;).oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOooo G4VPhysicalVolume* B4DetectorConstruction:DefineV olumes() / Geometry parametersG4int nofLayers = 10;定义有多少层 copynumberG4double absoThickness = 10.*mm;G4double gapThickness =5.*mm;G4double calorSizeXY = 10.*c

19、m;auto layerThickness = absoThickness + gapThickness;一层厚度auto calorThickness = nofLayers * layerThickness; 整体厚度auto worldSizeXY = 1.2 * calorSizeXY;auto worldSizeZ = 1.2 * calorThickness;/ Get materialsauto defaultMaterial = G4Material二GetMaterial(Galactic);auto absorberMaterial = G4Material:GetMate

20、rial(G4_Pb);auto gapMaterial = G4Material:GetMaterial(liquidArgon);if ( ! defaultMaterial | ! absorberMaterial | ! gapMaterial ) G4ExceptionDescription msg;msg Cannot retrieve materials already defined.;G4Exception(B4DetectorConstruction:DefineV olumes(), MyCode0001, FatalException, msg);)/ World/au

21、to worldS= new G4Box(World,/ its nameworldSizeXY/2, worldSizeXY/2, worldSizeZ/2); / its sizeauto worldLV=new G4LogicalV olume(worldS,/ its soliddefaultMaterial,/ its materialWorld);/ its nameauto worldPV=new G4PVPlacement(0,/ no rotationG4ThreeVector(),/ at (0,0,0)worldLV ,/ its logical volumeWorld,

22、/ its name0,/ its mother volumefalse,/ no boolean operation0,/ copy numberfCheckOverlaps);/ checking overlaps/ Calorimeter/auto calorimeters=new G4Box(Calorimeter, / its namecalorSizeXY/2, calorSizeXY/2, calorThickness/2); / its sizeauto calorLV=new G4LogicalV olume( calorimeterS, defaultMaterial, C

23、alorimeter);/ its solid/ its material/ 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/ Layer/auto layers=new G4Box(Layer,/ its namecalorSizeXY/2, calorSizeXY/2, layerThickness/2); / its

24、 sizeauto layerLV=new G4LogicalV olume( layerS,/ its soliddefaultMaterial, / its materialLayer);/ its namenew G4PVReplica(Layer,/its namelayerLV ,/its logical volumecalorLV ,/its motherkZAxis,/ axis of replicationnofLayers,/number of replicalayerThickness); / witdth of replica/ Absorber/auto absorbe

25、rS=new G4Box(Abso,/ its namecalorSizeXY/2, calorSizeXY/2, absoThickness/2); / its sizeauto absorberLV=new G4LogicalV olume(absorberS,/ its solidabsorberMaterial, / its material Abso);/ its namefAbsorberPV=new G4PVPlacement(0,/ no rotationG4ThreeVector(0., 0., -gapThickness/2), its position absorberL

26、V, / its logical volumeAbso, layerLV, false,0, fCheckOverlaps);Abso, layerLV, false,0, fCheckOverlaps);/ its mother volume/ no boolean operation/ copy number/ checking overlaps/ Gap/auto gapS=new G4Box(Gap,/ its namecalorSizeXY/2, calorSizeXY/2, gapThickness/2); / its sizeauto gapLV=new G4LogicalV o

27、lume( gapS,/ its solidgapMaterial, / its material Gap);/ its namefGapPV=new G4PVPlacement(0,/ no rotationG4ThreeVector(0., 0., absoThickness/2), / its position gapLV,/ its logical volumeGap,/ its namelayerLV ,/ its mother volumefalse,/ no boolean operation0,/ copy numberfCheckOverlaps); / checking o

28、verlaps/ print parameters /G4cout G4endl G4endl The calorimeter is nofLayers layers of: absoThickness/mm mm of GetName() + gapThickness/mm mm of GetName() G4endl SetVisAttributes (G4VisAttributes:GetInvisible();auto simpleBoxVisAtt= new G4VisAttributes(G4Colour(1.0,1.0,1.0);simpleBoxVisAtt-SetVisibi

29、lity(true);calorLV-SetVisAttributes(simpleBoxVisAtt);/ Always return the physical World /return worldPV; /oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOooovoid B4DetectorConstruction:ConstructSDandField() / Create global magnetic field messenger./ Uniform magnetic field is then created automatically if /

30、 the field value is not zero.G4ThreeVector fieldV alue;fMagFieldMessenger = new G4GlobalMagFieldMessenger(fieldValue); fMagFieldMessenger-SetVerboseLevel(1);/ Register the field messenger for deletingG4AutoDelete:Register(fMagFieldMessenger); 在 construct。函数中，定义了 Calorimeter、Layer、Abso以及Gap四个逻辑 体，它们之

31、间的关系是：构建一个大长方体 Calorimeter”,往里面复制并放置 10个规格相 同的长方体称之为层“ Layer，再往“ Layer”中分别放置“ Abso”以及“ Gap”。在 ConstructSDandField()函数体里，我们可以定义 SD探测器以及空间电磁场， 在本例 exampleB4a中我们只是定义了电磁场，可以通过run1.mac中命令行控制磁场属性：/globalField/setValue 0.2 0 0 tesla设置均匀磁场，强度 0.2特斯拉，方向沿 x轴exampleB4c中，我们在 ConstructSDandField()函数体里还定义了两个SD探测器

32、，这样做简化了后续数据抽取工作中对每一个step的位置判断，具体怎么操作到B4c例中再说。B4aActionInitialization.hh#ifndef B4aActionInitialization_h#define B4aActionInitialization_h 1#include G4VUserActionInitialization.hhclass B4DetectorConstruction;/ Action initialization class./class B4aActionInitialization : public G4VUserActionInitializa

33、tion public:B4aActionInitialization(B4DetectorConstruction*);virtual B4aActionInitialization();virtual void BuildForMaster() const;virtual void Build() const;private:B4DetectorConstruction* fDetConstruction; ;#endif在构造函数中，用到了 B4DetectorConstruction类的对象做形参，这是因为在后续判断 粒子当前位置过程中用到了探测器几何，用于获取几何的物理体，后面会说，

34、这里只需要知道用户初始化过程用到了B4DetectorConstruction类。B4aActionInitialization.ccB4aActionInitialization:B4aActionInitialization(B4DetectorConstruction* detConstruction): G4VUserActionInitialization(),fDetConstruction(detConstruction) ()/.oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOoooB4aActionInitialization二B4aActi

35、onInitialization() ()/.oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOooovoid B4aActionInitialization二BuildForMaster() const (SetUserAction(new B4RunAction); )/oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOooovoid B4aActionInitialization二Build() const (SetUserAction(new B4PrimaryGeneratorAction);SetUserAction(

36、new B4RunAction);auto eventAction = new B4aEventAction;SetUserAction(eventAction);SetUserAction(new B4aSteppingAction(fDetConstruction,eventAction); 在 step 里面用 到了 fDetConstruction 变量)将 detConstruction 赋值给 fDetConstruction，另外，对粒子发射器、RunAction、EventAction 进行用户行为设置，其中，B4aSteppingAction用到了探测器几何信息，一般是用于位

37、置判断。SetUserAction(new B4aSteppingAction(fDetConstruction,eventAction);B4PrimaryGeneratorAction.hhclass G4ParticleGun;class G4Event;/ The primary generator action class with particle gum./Ill It defines a single particle which hits the calorimeterIII perpendicular to the input face. The type of the p

38、articleIll can be changed via the G4 build-in commands of G4ParticleGun classIll (see the macros provided with this example).class B4PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction ( public:B4PrimaryGeneratorAction();virtual B4PrimaryGeneratorAction();virtual void GeneratePrimaries(G4E

39、vent* event);II set methodsvoid SetRandomFlag(G4bool value);private:G4ParticleGun* fParticleGun; II G4 particle gun);粒子发射器用的是 G4ParticleGun类。B4PrimaryGeneratorAction.cc#include B4PrimaryGeneratorAction.hh#include G4RunManager.hh#include G4LogicalV olumeStore.hh#include G4LogicalV olume.hh#include G4

40、Box.hh#include G4Event.hh#include G4ParticleGun.hh#include G4ParticleTable.hh#include G4ParticleDefinition.hh#include G4SystemOfUnits.hh#include Randomize.hh/.oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOoooB4PrimaryGeneratorAction:B4PrimaryGeneratorAction():G4VUserPrimaryGeneratorAction(), fParticleGun

41、(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.

42、,1.);fParticleGun-SetParticleEnergy(50.*MeV);/oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOoooB4PrimaryGeneratorAction:B4PrimaryGeneratorAction()delete fParticleGun;/oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOooovoid B4PrimaryGeneratorAction:GeneratePrimaries(G4Event* anEvent)/ This function is called at

43、the begining of event/ In order to avoid dependence of PrimaryGeneratorAction/ on DetectorConstruction class we get world volume/ from G4LogicalVolumeStore/G4double worldZHalfLength = 0.;auto worldLV = G4LogicalV olumeStore:GetInstance()-GetV olume(World); / 将G4LogicalVolumeStore 实例化，把存储在constructd的

44、逻辑体找出来/ Check that the world volume has box shapeG4Box* worldBox = nullptr;把 world solid 把 world solid 类worldBox = dynamic_cast(worldLV-GetSolid();型强制转化成box型if ( worldBox ) worldZHalfLength = worldBox-GetZHalfLength();抽取 word 的 Z 轴半长else G4ExceptionDescription msg;msg World volume of box shape not f

45、ound. G4endl;msg Perhaps you have changed geometry. G4endl;msg SetParticlePosition(G4ThreeVector(0., 0., -worldZHalfLength); 义在 word体上的 x-y面将源的位置定调用，开始fParticleGun-GeneratePrimaryVertex(anEvent);粒子枪的设置应该一样，设置一个默认初始值(可以在macro中更改)，默认是电子发射，电子能量50 MeV,方向沿z轴。在对源位置定义时，从world逻辑体中把实体 get出来后赋值给worldBox指针，此时w

46、orldBox指针的地址就是在几何体中构建的world实体，通过 worldZHalfLength = worldBox-GetZHalfLength() ;就是将 world 实体指向的内存空间的值 拿出来赋给 worldZHalfLength ,于是源的位置就出来了，是一个贴着world的x-y面的源。然后每一个 event 到来时调用 fParticleGun-GeneratePrimaryVertex(anEvent);发射粒子。B4RunAction.hhclass B4RunAction : public G4UserRunAction (public:B4RunAction();

47、virtual -B4RunAction();virtual void BeginOfRunAction(const G4Run*);virtual void EndOfRunAction(const G4Run*); ;B4RunAction继承于 G4UserRunAction基类，构造函数、析构函数以及每次run的begin和end虚方法。B4RunAction.cc (Begin)B4RunAction:B4RunAction(): G4UserRunAction()(/ set printing event number per each eventG4RunManager二GetR

48、unManager()-SetPrintProgress(1);/通过 RunManager 管理器获取进程，并打印/ Create analysis manager/ The choice of analysis technology is done via selectin of a namespace/ in B4Analysis.hhauto analysisManager = G4AnalysisManager:Instance(); / G4AnalysisManager 实例化对 象 analysisManagerG4cout Using GetType() SetHistoDi

49、rectoryName(histograms);/analysisManager-SetNtupleDirectoryName(ntuple);/analysisManager-SetVerboseLevel(1);analysisManager-SetNtupleMerging(true);/ Note: merging ntuples is available only with Root output/ Book histograms, ntuple/ Creating histogramsanalysisManager-CreateH1(Eabs,Edep in absorber, 1

50、00, 0., 800*MeV); 创建 histogramanalysisManager-CreateH1(Egap,Edep in gap, 100, 0., 100*MeV);analysisManager-CreateH1(Labs,trackL in absorber, 100, 0., 1*m); analysisManager-CreateH1(Lgap,trackL in gap, 100, 0., 50*cm);/ Creating ntuple创建 ntuple/ 分支/analysisManager-CreateNtuple(B4, Edep and TrackL); /

51、ntuple 名称： B4 analysisManager-CreateNtupleDColumn(Eabs);analysisManager-CreateNtupleDColumn(Egap); analysisManager-CreateNtupleDColumn(Labs); analysisManager-CreateNtupleDColumn(Lgap); analysisManager-FinishNtuple();/oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOoooB4RunAction:B4RunAction()delete G4Analy

52、sisManager:Instance();/oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOooovoid B4RunAction:BeginOfRunAction(const G4Run* /*run*/) /inform the runManager to save random number seed /G4RunManager:GetRunManager()-SetRandomNumberStore(true);/ Get analysis manager auto analysisManager = G4AnalysisManager:Instan

53、ce();/ Open an output file/G4String fileName = B4;analysisManager-OpenFile(fileName);打开B4打开B4文件在B4RunAction的构造函数中，通过runmanager获取进程，设置每隔多少event打印出 event数。analysisManager指针开始创建 histogram和ntuple用于数据存放， 相当于在root中 新建一个TH1F和tree。析构函数中，最后将删除analysisManager,释放内存空间。通过analysisManager打开B4”文件，用于在 EventAction中存储

54、信息。B4aEventAction.hhclass B4aEventAction : public G4UserEventActionpublic:B4aEventAction();virtual B4aEventAction();virtual void BeginOfEventAction(const G4Event* event);virtual void EndOfEventAction(const G4Event* event);void AddAbs(G4double de, G4double dl);/函数，吸收体属性void AddGap(G4double de, G4doub

55、le dl);/函数，间隙属性private:G4double fEnergyAbs;G4double fEnergyGap;G4double fTrackLAbs;G4double fTrackLGap;/ inline functionsinline void B4aEventAction:AddAbs(G4double de, G4double dl) fEnergyAbs += de;fTrackLAbs += dl;inline void B4aEventAction:AddGap(G4double de, G4double dl) fEnergyGap += de;fTrackLG

56、ap += dl;)/.oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOooo#endif在 B4aEventAction 类中自定义了两个内联系函数，AddAbs(G4double de, G4double dl)和AddGap(G4double de, G4double dl),在类体外对他们进行了定义，分别用于累加。B4aEventAction.cc (Begin)B4aEventAction:B4aEventAction(): G4UserEventAction(),fEnergyAbs(0.),fEnergyGap(0.),fTrackLAbs(0.

57、),fTrackLGap(0.) )/oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOoooB4aEventAction:B4aEventAction() )/oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOooovoid B4aEventAction:BeginOfEventAction(const G4Event* /*event*/) / initialisation per event/每一个 event 开始，首先初始化fEnergyAbs = 0.;fEnergyGap = 0.;fTrackLAbs = 0.;f

58、TrackLGap = 0.; )构造函数中，对私有变量 fEnergyAbs、fEnergyGap、fTrackLAbs、fTrackLGap 赋初值 0, 在BeginOfEvebtAction()中，每一次 event开始，都对他们进行初始化。B4aSteppingAction.hh#ifndef B4aSteppingAction_h#define B4aSteppingAction_h 1#include G4UserSteppingAction.hhclass B4DetectorConstruction;class B4aEventAction;/ Stepping action

59、 class./ In UserSteppingAction() there are collected the energy deposit and track/ lengths of charged particles in Absober and Gap layers and/ updated in B4aEventAction.class B4aSteppingAction : public G4UserSteppingActionpublic:B4aSteppingAction(const B4DetectorConstruction* detectorConstruction,与探

60、测器的位置是有关系的B4aEventAction* eventAction);virtual B4aSteppingAction();virtual void UserSteppingAction(const G4Step* step);private:const B4DetectorConstruction* fDetConstruction;B4aEventAction* fEventAction;.oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOooo#endif在 B4aSteppingAction 的构造函数中， B4aSteppingAction(