蒙特卡洛模擬GEANT4例子exampleB4a源碼解讀

1、Geant4 ExampleB4a 源碼解讀今天繼續(xù)簡單講一下exampleB4中的幾個子例子，先是exampleB4a，仍然是采用Steppi ngAction 步一步抽取事件，這種方法的好處是能夠很詳細的考慮、輸出每一步的信 息，但有時候我們只希望對感興趣部分（位置）進行事件抽取，則用SD和Hit要明顯方便一些，例如在exampleB4c中用到的就是 SD和Hit進行數(shù)據(jù)抽取。exampleB4a和exampleB4c 描述的是同一個幾何、物理過程，只是數(shù)據(jù)抽取的方式不同，因此將這兩個例子放在一起具 有比較意義。本次先對 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&

3、quot;#en dif#in clude "G4UIma nager.hh"#in clude "G4UIcomma nd.hh"#include "FTFP_BERT.hh" #include "Randomize.hh"#include "G4VisExecutive.hh"#include "G4UIExecutive.hh"/oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOooonamespace void PrintUsage()

4、 G4cerr << " Usage: " << G4endl;G4cerr << " exampleB4a -m macro -u UIsession -t nThreads" << G4endl;G4cerr << " note: -t option is available only for multi-threaded mode." << G4endl;/oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOooo int ma

5、in(int argc,char* argv) / Evaluate arguments/if ( argc > 7 ) PrintUsage(); return 1;G4String macro;G4String session;#ifdef G4MULTITHREADEDG4int nThreads = 0;#endiffor ( G4int i=1; i<argc; i=i+2 ) if( G4String(argvi) = "-m" ) macro = argvi+1;else if ( G4String(argvi) = "-u"

6、) session = argvi+1;#ifdef G4MULTITHREADEDelse if ( G4String(argvi) = "-t" ) nThreads = G4UIcommand:ConvertToInt(argvi+1);#endif else PrintUsage();return 1;/ Detect interactive mode (if no macro provided) and define UI session/G4UIExecutive* ui = nullptr;if ( ! macro.size() ) ui = new G4UI

7、Executive(argc, argv, session);/ Optionally: choose a different Random engine./ G4Random:setTheEngine(new CLHEP:MTwistEngine);/ Construct the default run manager/#ifdef G4MULTITHREADEDauto runManager = new G4MTRunManager;if ( nThreads > 0 ) runManager->SetNumberOfThreads(nThreads);#elseauto ru

8、nManager = new G4RunManager;#endif/ Set mandatory initialization classes/auto detConstruction = new B4DetectorConstruction(); runManager->SetUserInitialization(detConstruction);auto physicsList = new FTFP_BERT; runManager->SetUserInitialization(physicsList);auto actionInitialization = new B4aA

9、ctionInitialization(detConstruction); runManager->SetUserInitialization(actionInitialization);/ Initialize visualization/auto visManager = new G4VisExecutive;/ G4VisExecutive can take a verbosity argument - see /vis/verbose guidance./ G4VisManager* visManager = new G4VisExecutive("Quiet"

10、;); visManager->Initialize();/ Get the pointer to the User Interface managerauto UImanager = G4UImanager:GetUIpointer();/ Process macro or start UI session/if ( macro.size() ) / batch modeG4String command = "/control/execute "UImanager->ApplyCommand(command+macro);else / interactive

11、mode : define UI sessionUImanager->ApplyCommand("/control/execute init_vis.mac");if (ui->IsGUI() UImanager->ApplyCommand("/control/execute gui.mac");ui->SessionStart();delete ui;/ Job termination/ Free the store: user actions, physics_list and detector_description are

12、/ owned and deleted by the run manager, so they should not be deleted/ in the main() program !delete visManager;delete runManager;用戶在終端輸入字符串個數(shù)大于 7 則調(diào)用 PrintUsage 提示出錯，并返回 1，跳出 main()函數(shù)。通過for循環(huán)對輸入字符串掃描識別，若程序名后字符串為"-m”給macro賦值該字符串后一個字符串；若程序名后字符串為“ -u”給session賦值該字符串后一個字符串，否則調(diào)用PrintUsage提示出錯，并返回 1跳

13、出main()函數(shù)。定義一個 UI操作類對象，ui指針 初始化空指針，若用戶沒有調(diào)用 macro 控制文件，則 new 一個 ui 指針，分內(nèi)存空間辟。初 始化探測器、物理過程以及用戶行為，其中，用戶行為與探測器構(gòu)造有關(guān)。new 一個可視化管理器并對其初始化。定義一個 UI 管理類，通過判斷條件 macro 是否為空，來決定執(zhí)行交 互式界面模式還是批處理模式。最后delete所有管理類，釋放內(nèi)存。B4DetectorConstruction.hh#ifndef B4DetectorConstruction_h #define B4DetectorConstruction_h 1#include

14、 "G4VUserDetectorConstruction.hh"#include "globals.hh"class G4VPhysicalVolume;class G4GlobalMagFieldMessenger;class B4DetectorConstruction : public G4VUserDetectorConstructionpublic:B4DetectorConstruction();virtual B4DetectorConstruction();public:virtual G4VPhysicalV olume* Const

15、ruct();virtual void ConstructSDandField();/ get methods/const G4VPhysicalVolume* GetAbsorberPV() const; /const 修飾指針， 指針的內(nèi)存空間 數(shù)據(jù)不能改變；函數(shù)重載const G4VPhysicalV olume* GetGapPV() const;private:/ methods/void DefineMaterials();G4VPhysicalV olume* DefineV olumes();/ data members/static G4ThreadLocal G4Globa

16、lMagFieldMessenger* fMagFieldMessenger;/ magnetic field messengerG4VPhysicalV olume* fAbsorberPV; / the absorber physical volumeG4VPhysicalV olume* fGapPV; / the gap physical volumeG4bool fCheckOverlaps; / option to activate checking of volumes overlaps;/ inline functionsinline const G4VPhysicalV ol

17、ume* B4DetectorConstruction:GetAbsorberPV() const return fAbsorberPV;inline const G4VPhysicalV olume* B4DetectorConstruction:GetGapPV() const return fGapPV;B4DetectorConstruction 繼承于 G4VUserDetectorConstruction 基類， 聲明構(gòu)造函數(shù)、 析構(gòu) 函數(shù)，Con struct。返回物理體。虛方法Con structSDa ndField()用于定義SD探測器以及空間內(nèi) 的電磁場。 GetAbso

18、rberPV() 和 GetGapPV() 為自定義的內(nèi)聯(lián)函數(shù)在類體外進行定義，分別返 回物理體 fAbsorberPV 和 fGapPV (私有數(shù)據(jù)成員) 。 G4GlobalMagFieldMessenger 類用于描 述電磁場的大小、方向等信息。B4DetectorConstruction.ccB4DetectorConstruction:B4DetectorConstruction(): G4VUserDetectorConstruction(),fAbsorberPV(nullptr),fGapPV(nullptr),fCheckOverlaps(true)/oooOO0OOoooo

19、ooOO0OOooooooOO0OOooooooOO0OOoooB4DetectorConstruction:B4DetectorConstruction()/oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOoooG4VPhysicalVolume* B4DetectorConstruction:Construct() / Define materialsDefineMaterials();/ Define volumesreturn DefineVolumes();/oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOooovo

20、id 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 protons;G4double density;new G4M

21、aterial("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 materialsG4cout

22、<< *(G4Material:GetMaterialTable() << G4endl;/oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOoooG4VPhysicalVolume* B4DetectorConstruction:DefineV olumes()/ Geometry parametersG4int nofLayers = 10;/ 定義有多少層 copynumberG4double absoThickness = 10.*mm;G4double gapThickness = 5.*mm;G4double calorSiz

23、eXY = 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("Galactic");auto absorberMater

24、ial = G4Material:GetMaterial("G4_Pb");auto gapMaterial = G4Material:GetMaterial("liquidArgon");if ( ! defaultMaterial | ! absorberMaterial | ! gapMaterial ) G4ExceptionDescription msg;msg << "Cannot retrieve materials already defined."G4Exception("B4DetectorC

25、onstruction:DefineV olumes()", "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 G4PVPlac

26、ement(0,G4ThreeVector(), worldLV , "World", 0, false, 0, fCheckOverlaps);/ 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(

27、"Calorimeter",calorSizeXY/2, calorSizeXY/2, calorThickness/2); / its sizeauto calorLV= new G4LogicalV olume(calorimeterS, defaultMaterial, "Calorimeter");new G4PVPlacement(0,G4ThreeVector(), calorLV, "Calorimeter", worldLV , false, 0, fCheckOverlaps);/ its solid/ its ma

28、terial/ 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 sizeauto layerLV= new G4LogicalV olume(la

29、yerS,/ its soliddefaultMaterial, / its material "Layer");/ its namenew G4PVReplica("Layer",/ its namelayerLV,/ its logical volumecalorLV,/ its motherkZAxis,/ axis of replicationnofLayers,/ number of replicalayerThickness); / witdth of replica/ Absorber/auto absorberS= new G4Box(&

30、quot;Abso",/ its namecalorSizeXY/2, calorSizeXY/2, absoThickness/2); / its sizeauto absorberLV= new G4LogicalV olume(absorberS, / its solid absorberMaterial, / its material "Abso");/ its namefAbsorberPV= new G4PVPlacement(0,/ no rotationG4ThreeVector(0., 0., -gapThickness/2), / its po

31、sition absorberLV,/ its logical volume"Abso", layerLV, false,0, fCheckOverlaps);/ its name/ 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/ its

32、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"Gap",/ its namelayerLV,/ its mother volumefalse,/ no boolean operation0,/ copy

33、 numberfCheckOverlaps); / checking overlaps/ print parameters/G4cout<< G4endl<< "" << G4endl<< "-> The calorimeter is " << nofLayers << " layers of: "<< absoThickness/mm << "mm of " << absorberMaterial-&

34、gt;GetName() << " + "<< gapThickness/mm << "mm of " << gapMaterial->GetName() << " " << G4endl<< "" << G4endl;/ Visualization attributes/ worldLV->SetVisAttributes (G4VisAttributes:GetInvisible();auto simpl

35、eBoxVisAtt= new G4VisAttributes(G4Colour(1.0,1.0,1.0); simpleBoxVisAtt->SetVisibility(true);calorLV->SetVisAttributes(simpleBoxVisAtt);/ Always return the physical World/return worldPV;/oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOooovoid B4DetectorConstruction:ConstructSDandField()/ Create global

36、 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 deleting G4AutoDelet

37、e:Register(fMagFieldMessenger);在 con struct。函數(shù)中，定義了" Calorimeter ”、"Layer ”、“ Abso ” 以及"Gap ”四個邏輯 體，它們之間的關(guān)系是：構(gòu)建一個大長方體"Calorimeter ”，往里面復制并放置 10個規(guī)格相同的長方體稱之為層"Layer ”，再往"Layer ”中分別放置" Abso ”以及"Gap ”。在ConstructSDandField() 函數(shù)體里，我們可以定義 SD 探測器以及空間電磁場， 在本例 exampleB4

38、a中我們只是定義了電磁場，可以通過run 1.mac中命令行控制磁場屬性：/globalField/setValue 0.2 0 0 tesla/設(shè)置均勻磁場，強度0.2特斯拉，方向沿 x軸exampleB4c 中，我們在 ConstructSDandField() 函數(shù)體里還定義了兩個 SD 探測器，這樣 做簡化了后續(xù)數(shù)據(jù)抽取工作中對每一個 step 的位置判斷，具體怎么操作到 B4c 例中再說。B4aActionInitialization.hh#ifndef B4aActionInitialization_h#define B4aActionInitialization_h 1#incl

39、ude "G4VUserActionInitialization.hh"class B4DetectorConstruction;/ Action initialization class./class B4aActionInitialization : public G4VUserActionInitializationpublic:B4aActionInitialization(B4DetectorConstruction*);virtual B4aActionInitialization();virtual void BuildForMaster() const;vi

40、rtual void Build() const;private:B4DetectorConstruction* fDetConstruction;#endif在構(gòu)造函數(shù)中，用到了 B4DetectorConstruction 類的對象做形參，這是因為在后續(xù)判斷 粒子當前位置過程中用到了探測器幾何， 用于獲取幾何的物理體， 后面會說， 這里只需要知 道用戶初始化過程用到了 B4DetectorConstruction 類。B4aActionInitialization.ccB4aActionInitialization:B4aActionInitialization(B4DetectorCon

41、struction* detConstruction): G4VUserActionInitialization(),fDetConstruction(detConstruction) /oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOoooB4aActionInitialization:B4aActionInitialization() /oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOooovoid B4aActionInitialization:BuildForMaster() constSetUserAction(ne

42、w B4RunAction);/oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOooovoid B4aActionInitialization:Build() const SetUserAction(new B4PrimaryGeneratorAction);SetUserAction(new B4RunAction);auto eventAction = new B4aEventAction;SetUserAction(eventAction);SetUserAction(new B4aSteppingAction(fDetConstruction,even

43、tAction);/在 step 里面用到了 fDetConstruction 變量將 detConstruction 賦值給 fDetConstruction ，另外，對粒子發(fā)射器、 RunAction 、EventAction 進行用戶行為設(shè)置，其中， B4aSteppingAction 用到了探測器幾何信息，一般是用于位置判 斷。SetUserAction(new B4aSteppingAction(fDetConstruction,eventAction);B4PrimaryGeneratorAction.hhclass G4ParticleGun;class G4Event;/ Th

44、e primary generator action class with particle 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 B4PrimaryGen

45、eratorAction : public G4VUserPrimaryGeneratorAction public:B4PrimaryGeneratorAction();virtual B4PrimaryGeneratorAction();virtual void GeneratePrimaries(G4Event* event);/ set methodsvoid SetRandomFlag(G4bool value);private:G4ParticleGun* fParticleGun; / G4 particle gun;粒子發(fā)射器用的是 G4ParticleGun 類。B4Prim

46、aryGeneratorAction.cc#include "B4PrimaryGeneratorAction.hh"#include "G4RunManager.hh"#include "G4LogicalV olumeStore.hh"#include "G4LogicalV olume.hh"#include "G4Box.hh"#include "G4Event.hh"#include "G4ParticleGun.hh"#include &quo

47、t;G4ParticleTable.hh"#include "G4ParticleDefinition.hh"#include "G4SystemOfUnits.hh"#include "Randomize.hh"/oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOoooB4PrimaryGeneratorAction:B4PrimaryGeneratorAction(): G4VUserPrimaryGeneratorAction(),fParticleGun(nullptr)G4int n

48、ofParticles = 1;fParticleGun = new G4ParticleGun(nofParticles);/ default particle kinematic/auto particleDefinition= G4ParticleTable:GetParticleTable()->FindParticle("e-"); fParticleGun->SetParticleDefinition(particleDefinition); fParticleGun->SetParticleMomentumDirection(G4ThreeV

49、ector(0.,0.,1.); fParticleGun->SetParticleEnergy(50.*MeV);/oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOoooB4PrimaryGeneratorAction:B4PrimaryGeneratorAction()delete fParticleGun;/oooOO0OOooooooOO0OOooooooOO0OOooooooOO0OOooovoid B4PrimaryGeneratorAction:GeneratePrimaries(G4Event* anEvent)/ This functi

50、on is called at 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"); G4Logica

51、lVolumeStore 實例化，把存儲在 constructd 的邏輯體找出來/ Check that the world volume has box shapeG4Box* worldBox = nullptr;if ( worldLV ) worldBox = dynamic_cast<G4Box*>(worldLV->GetSolid(); 型 強制轉(zhuǎn)化成 box 型/把 world solid 類if ( worldBox ) worldZHalfLength = worldBox->GetZHalfLength();/ 抽取 word 的 Z 軸半長els

52、e G4ExceptionDescription msg;msg << "World volume of box shape not found." << G4endl; msg << "Perhaps you have changed geometry." << G4endl; msg << "The gun will be place in the center."G4Exception("B4PrimaryGeneratorAction:GeneratePri

53、maries()", "MyCode0002", JustWarning, msg);/ Set gun position fParticleGun->SetParticlePosition(G4ThreeVector(0., 0., -worldZHalfLength);/將源的位置定義在 word 體上的 x-y 面fParticleGun->GeneratePrimaryVertex(anEvent);/調(diào)用，開始粒子槍的設(shè)置應該一樣，設(shè)置一個默認初始值(可以在 macro 中更改)，默認是電子發(fā) 射，電子能量 50 MeV ，方向沿 z 軸。在

54、對源位置定義時， 從 world 邏輯體中把實體 get 出來 后賦值給 worldBox 指針， 此時 worldBox 指針的地址就是在幾何體中構(gòu)建的 world 實體， 通 過 worldZHalfLength = worldBox->GetZHalfLength() ；就是將 world 實體指向的內(nèi)存空間的值 拿出來賦給 worldZHalfLength ，于是源的位置就出來了，是一個貼著 world 的 x-y 面的源。然后每一個 event 到來時調(diào)用 fParticleGun->GeneratePrimaryVertex(anEvent) ；發(fā)射粒子。B4RunAc

55、tion.hhclass B4RunAction : public G4UserRunActionpublic:B4RunAction();virtual B4RunAction();virtual void BeginOfRunAction(const G4Run*);virtual voidEndOfRunAction(const G4Run*);B4RunAction 繼承于 G4UserRunAction 基類，構(gòu)造函數(shù)、 析構(gòu)函數(shù)以及每次 run 的 begin 和 end 虛方法。B4RunAction.cc (Begin)B4RunAction:B4RunAction(): G4

56、UserRunAction() / set printing event number per each eventG4RunManager:GetRunManager()->SetPrintProgress(1);/通過 RunManager 管理器獲取進程，并打印/ Create analysis manager/ The choice of analysis technology is done via selectin of a namespace/ in B4Analysis.hhauto analysisManager = G4AnalysisManager:Instance

57、(); / G4AnalysisManager 實例化對 象 analysisManagerG4cout << "Using " << analysisManager->GetType() << G4endl;/ Create directories /analysisManager->SetHistoDirectoryName("histograms");/analysisManager->SetNtupleDirectoryName("ntuple"); / analysisM

58、anager->SetVerboseLevel(1);analysisManager->SetNtupleMerging(true);/ Note: merging ntuples is available only with Root output/ Book histograms, ntuple/ Creating histogramsanalysisManager->CreateH1("Eabs","Edep in absorber", 100, 0., 800*MeV);/ 創(chuàng) 建histogramanalysisManager->CreateH1("Egap","Edep in gap", 100, 0., 100*MeV); analysisManager->CreateH1("Labs