geant4如何简单地获取粒子在闪烁体上产生的脉冲波形,请给我完整的代码
时间: 2024-02-01 09:14:20 浏览: 202
获取粒子在闪烁体上产生的脉冲波形需要进行光学模拟,需要使用Geant4中的光学模拟工具。以下是一个获取闪烁体上光子脉冲波形的示例代码:
#include "G4RunManager.hh"
#include "G4UImanager.hh"
#include "G4UIterminal.hh"
#include "G4UItcsh.hh"
#include "G4VisExecutive.hh"
#include "G4UIExecutive.hh"
#include "G4LogicalVolume.hh"
#include "G4Material.hh"
#include "G4NistManager.hh"
#include "G4Box.hh"
#include "G4SDManager.hh"
#include "G4VSensitiveDetector.hh"
#include "G4OpticalPhoton.hh"
#include "G4Scintillation.hh"
#include "G4OpBoundaryProcess.hh"
#include "G4LogicalBorderSurface.hh"
#include "G4LogicalSkinSurface.hh"
#include "G4PhysicalConstants.hh"
#include "G4SystemOfUnits.hh"
#include "G4ios.hh"
#include <fstream>
using namespace std;
class DetectorSD : public G4VSensitiveDetector
{
public:
DetectorSD(G4String name) : G4VSensitiveDetector(name) {}
virtual ~DetectorSD() {}
virtual void Initialize(G4HCofThisEvent* /*HCE*/)
{
fPhotonCount = 0;
}
virtual G4bool ProcessHits(G4Step* aStep, G4TouchableHistory* /*ROhist*/)
{
G4Track* aTrack = aStep->GetTrack();
if (aTrack->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition())
{
fPhotonCount++;
// 记录光子入射闪烁体的时间
fPhotonTimes.push_back(aTrack->GetGlobalTime());
}
return true;
}
vector<G4double> GetPhotonTimes()
{
return fPhotonTimes;
}
private:
G4int fPhotonCount;
vector<G4double> fPhotonTimes;
};
int main(int argc, char** argv)
{
G4RunManager* runManager = new G4RunManager;
// 初始化随机数生成器
G4Random::setTheEngine(new CLHEP::RanecuEngine);
G4long seed = time(NULL);
G4Random::setTheSeed(seed);
// 定义材料
G4NistManager* nistManager = G4NistManager::Instance();
G4Material* matScint = nistManager->FindOrBuildMaterial("G4_PLASTIC_SC_VINYLTOLUENE");
// 定义实验室空气
G4double labAirDensity = 1.29e-3 * g/cm3;
G4Material* matAir = new G4Material("Air", 1.29e-3 * g/cm3, 2);
matAir->AddElement(nistManager->FindOrBuildElement("N"), 70 * perCent);
matAir->AddElement(nistManager->FindOrBuildElement("O"), 30 * perCent);
// 定义闪烁体
G4double scintSizeXY = 10 * cm;
G4double scintSizeZ = 1 * cm;
G4Box* solidScint = new G4Box("scint", scintSizeXY/2, scintSizeXY/2, scintSizeZ/2);
G4LogicalVolume* logicScint = new G4LogicalVolume(solidScint, matScint, "scint");
// 定义光学性质
G4double photonEnergy[] = { 2.034*eV, 4.136*eV };
G4double rIndexScint[] = { 1.58, 1.58 };
G4double absorptionLength[] = { 4.16*m, 4.16*m };
G4double scintilFast[] = { 1.00, 0.00 };
G4double scintilSlow[] = { 0.00, 1.00 };
G4MaterialPropertiesTable* mptScint = new G4MaterialPropertiesTable();
mptScint->AddProperty("RINDEX", photonEnergy, rIndexScint, 2);
mptScint->AddProperty("ABSLENGTH", photonEnergy, absorptionLength, 2);
mptScint->AddProperty("FASTCOMPONENT", photonEnergy, scintilFast, 2);
mptScint->AddProperty("SLOWCOMPONENT", photonEnergy, scintilSlow, 2);
mptScint->AddConstProperty("SCINTILLATIONYIELD", 10000./MeV);
mptScint->AddConstProperty("RESOLUTIONSCALE", 1.0);
mptScint->AddConstProperty("FASTTIMECONSTANT", 1.*ns);
mptScint->AddConstProperty("SLOWTIMECONSTANT", 10.*ns);
mptScint->AddConstProperty("YIELDRATIO", 1.0);
logicScint->SetMaterialPropertiesTable(mptScint);
// 定义探测器敏感区域
DetectorSD* detectorSD = new DetectorSD("Detector");
G4SDManager* sdManager = G4SDManager::GetSDMpointer();
sdManager->AddNewDetector(detectorSD);
logicScint->SetSensitiveDetector(detectorSD);
// 定义闪烁体表面
G4OpticalSurface* opScintSurface = new G4OpticalSurface("ScintSurface");
opScintSurface->SetType(dielectric_dielectric);
opScintSurface->SetFinish(polished);
opScintSurface->SetModel(unified);
G4LogicalSkinSurface* scintSurface = new G4LogicalSkinSurface(
"ScintSurface", logicScint, opScintSurface);
// 定义实验室
G4double worldSizeXY = 20 * cm;
G4double worldSizeZ = 20 * cm;
G4Box* solidWorld = new G4Box("world", worldSizeXY/2, worldSizeXY/2, worldSizeZ/2);
G4LogicalVolume* logicWorld = new G4LogicalVolume(solidWorld, matAir, "world");
// 定义实验室边界
G4double borderSize = 1 * mm;
G4Box* solidBorder = new G4Box("border", worldSizeXY/2+borderSize, worldSizeXY/2+borderSize, worldSizeZ/2+borderSize);
G4LogicalVolume* logicBorder = new G4LogicalVolume(solidBorder, matAir, "border");
new G4PVPlacement(0, G4ThreeVector(), logicBorder, "border", logicWorld, false, 0);
new G4LogicalBorderSurface("worldToBorder", logicWorld, logicBorder, opScintSurface);
// 定义实验室中的光源
G4double tubeRadius = 0.5 * cm;
G4double tubeLength = 10 * cm;
G4Tubs* solidTube = new G4Tubs("tube", 0, tubeRadius, tubeLength/2, 0, 360*deg);
G4LogicalVolume* logicTube = new G4LogicalVolume(solidTube, matAir, "tube");
new G4PVPlacement(0, G4ThreeVector(0, 0, -5*cm-tubeLength/2), logicTube, "tube", logicWorld, false, 0);
// 定义光学性质
G4double rIndexAir[] = { 1.0, 1.0 };
G4double absorptionLengthAir[] = { 1e10*m, 1e10*m };
G4MaterialPropertiesTable* mptAir = new G4MaterialPropertiesTable();
mptAir->AddProperty("RINDEX", photonEnergy, rIndexAir, 2);
mptAir->AddProperty("ABSLENGTH", photonEnergy, absorptionLengthAir, 2);
G4OpticalSurface* opAirSurface = new G4OpticalSurface("AirSurface");
opAirSurface->SetType(dielectric_dielectric);
opAirSurface->SetFinish(polished);
opAirSurface->SetModel(unified);
G4LogicalSkinSurface* airSurface = new G4LogicalSkinSurface(
"AirSurface", logicTube, opAirSurface);
opAirSurface->SetMaterialPropertiesTable(mptAir);
// 定义光学过程
G4ProcessManager* pm = G4OpticalPhoton::OpticalPhoton()->GetProcessManager();
pm->AddProcess(new G4Scintillation(), 0, 1, 1);
pm->AddProcess(new G4OpBoundaryProcess(), 0, -1, 2);
// 定义可视化
G4VisManager* visManager = new G4VisExecutive;
visManager->Initialize();
// 开始模拟
runManager->SetUserInitialization(logicWorld);
runManager->Initialize();
G4UImanager* uiManager = G4UImanager::GetUIpointer();
G4UIExecutive* ui = new G4UIExecutive(argc, argv);
uiManager->ApplyCommand("/control/execute vis.mac");
ui->SessionStart();
// 获取光子入射闪烁体的时间
vector<G4double> photonTimes = detectorSD->GetPhotonTimes();
sort(photonTimes.begin(), photonTimes.end());
// 生成脉冲波形
G4double timeStep = 0.1 * ns;
G4double pulseWidth = 20 * ns;
G4double pulseAmplitude = 1000;
G4double pulseBaseline = 100;
ofstream outFile("pulse.dat");
for (G4double t = photonTimes[0]; t < photonTimes.back(); t += timeStep)
{
G4double pulse = pulseBaseline;
for (G4double ti : photonTimes)
{
if (ti >= t && ti < t+pulseWidth)
{
pulse += pulseAmplitude * exp(-(ti-t)/pulseWidth);
}
}
outFile << t/ns << " " << pulse << endl;
}
outFile.close();
// 清理内存
delete visManager;
delete ui;
delete runManager;
delete matAir;
delete opScintSurface;
delete opAirSurface;
delete nistManager;
return 0;
}
在这个示例中,我们定义了一个 DetectorSD
类作为探测器敏感区域。当光子进入该区域时,会将光子的时间记录下来。在模拟结束后,我们可以获取所有光子在闪烁体上产生的时间,并根据这些时间生成脉冲波形。具体实现细节可以参考代码中的注释。
相关推荐













