Geant4 Version: geant4-v11.3.0
Operating System: Windows 11 64bit
Compiler/Version: Visual Studio 2022
CMake Version: 4.0.1
Hello everyone,
I try to simulate the response of scintillation detector (EJ276) using Cs-137 as gamma source and when i run the simulation i see that few interaction of particles via Compton and photoelectric ( In PMT). I want to know the issue. Is it related to geometry or Physicslist.
This is the Geometry i have.
void DetectorConstruction::SetupGeometry() {
// World volume - large enough to contain everything
const G4double worldSize = 3.0 * (fScSize + fPMTLength + fLightGuideLength + 10.0 * cm);
auto* sWorld = new G4Box("World", worldSize / 2, worldSize / 2, worldSize / 2);
fWorldLV = new G4LogicalVolume(sWorld, fWorldMat, "WorldLV");
fWorldLV->SetVisAttributes(fVisAttributes[0]);
fWorldPV = new G4PVPlacement(nullptr, {}, fWorldLV, "WorldPV", nullptr, false, 0);
// 1. Dark box (light-tight enclosure) - placed in world
const G4double darkBoxThickness = 2.0 * mm;
const G4double darkBoxSizeXY = fScSize + 2 * 0.5 * mm + 2 * darkBoxThickness;
const G4double darkBoxSizeZ = fScSize + (fPMTLength + fLightGuideLength + 0.1 * mm + darkBoxThickness);
auto* sDarkBox = new G4Box("DarkBox", darkBoxSizeXY / 2, darkBoxSizeXY / 2, darkBoxSizeZ / 2);
G4LogicalVolume* darkBoxLV = new G4LogicalVolume(sDarkBox, fAlFoil, "DarkBoxLV");
darkBoxLV->SetVisAttributes(G4VisAttributes(G4Colour(0.1, 0.1, 0.1))); // Dark gray
// Place dark box in world at origin
G4VPhysicalVolume* darkBoxPV = new G4PVPlacement(nullptr, G4ThreeVector(0, 0, 0), darkBoxLV, "DarkBoxPV", fWorldLV, false, 0);
// 2. Wrapping layer (0.5 mm thick) - placed in dark box
const G4double wrapThickness = 0.5 * mm;
// Wrap should be large enough to contain scintillator and front assembly
const G4double wrapZ = fScSize + (fPMTLength + fLightGuideLength + 0.1 * mm);
auto* sWrap = new G4Box("Wrap",
fScSize / 2 + wrapThickness,
fScSize / 2 + wrapThickness,
wrapZ / 2);
G4Material* wrapMat = fWrapType == "Teflon" ? fTeflon :
fWrapType == "ESR" ? fESR :
fWrapType == "AlFoil" ? fAlFoil : fTiO2Paint;
fWrapLV = new G4LogicalVolume(sWrap, wrapMat, "WrapLV");
fWrapLV->SetVisAttributes(fVisAttributes[2]);
// Place wrap in dark box - offset to have scintillator at one end
fWrapPV = new G4PVPlacement(nullptr, G4ThreeVector(0, 0, -fScSize / 2), fWrapLV, "WrapPV", darkBoxLV, false, 0);
// 3. Scintillator - placed at one end of wrap
auto* sScint = new G4Box("Scintillator", fScSize / 2, fScSize / 2, fScSize / 2);
fScintillatorLV = new G4LogicalVolume(sScint, fScintMat, "ScintillatorLV");
fScintillatorLV->SetVisAttributes(fVisAttributes[1]);
// Place scintillator at the negative Z end of wrap (one side only)
fScintillatorPV = new G4PVPlacement(nullptr,
G4ThreeVector(0, 0, -wrapZ / 2 + fScSize / 2), // At negative end of wrap
fScintillatorLV, "ScintillatorPV", fWrapLV, false, 0);
// 4. Grease layer (0.1 mm)
const G4double greaseThickness = 0.1 * mm;
auto* sGrease = new G4Tubs("GreaseDisk", 0, fPMTRadius, greaseThickness / 2, 0, twopi);
fGreaseLV = new G4LogicalVolume(sGrease, fGrease, "GreaseLV");
fGreaseLV->SetVisAttributes(fVisAttributes[5]);
// 5. Light guide (if enabled)
if (fUseLightGuide) {
// Conical light guide to match scintillator face to PMT diameter
const G4double r1 = fScSize / 2; // Match scintillator face
const G4double r2 = fPMTRadius; // Match PMT diameter
auto* sLightGuide = new G4Cons("LightGuide",
0, r1, 0, r2, fLightGuideLength / 2, 0, twopi);
fLightGuideLV = new G4LogicalVolume(sLightGuide, fPMMA, "LightGuideLV");
fLightGuideLV->SetVisAttributes(fVisAttributes[4]);
}
// 6. PMT assembly
auto* sPMT = new G4Tubs("PMT_tube", 0, fPMTRadius, fPMTLength / 2, 0, twopi);
fPMTLV = new G4LogicalVolume(sPMT, fPMTGlass, "PMTLV");
fPMTLV->SetVisAttributes(fVisAttributes[3]);
// 7. Photocathode (thin layer inside PMT)
const G4double photoThickness = 0.05 * mm;
const G4double photoRadius = fPMTRadius - 0.5 * mm; // Slightly smaller than PMT
auto* sPhoto = new G4Tubs("PhotocathodeDisk", 0, photoRadius, photoThickness / 2, 0, twopi);
fPhotoLV = new G4LogicalVolume(sPhoto, fBialkali, "PhotocathodeLV");
fPhotoLV->SetVisAttributes(G4VisAttributes(G4Colour(0.0, 0.0, 1.0, 0.7))); // Blue
// 8. PMT magnetic shield (mu-metal)
const G4double shieldThickness = 1.0 * mm;
auto* sShield = new G4Tubs("PMTShield",
fPMTRadius,
fPMTRadius + shieldThickness,
fPMTLength / 2 + shieldThickness,
0, twopi);
G4LogicalVolume* shieldLV = new G4LogicalVolume(sShield, fMuMetal, "ShieldLV");
shieldLV->SetVisAttributes(G4VisAttributes(G4Colour(0.5, 0.3, 0.1))); // Bronze
// Positioning INSIDE WRAP - only on positive Z side of scintillator
// Starting from the center of the wrap volume
G4double zPos = -wrapZ / 2 + fScSize / 2; // Start at negative end of scintillator
// Grease layer - placed at positive end of scintillator
G4double zGrease = zPos + fScSize / 2 + greaseThickness / 2;
fGreaseFrontPV = new G4PVPlacement(nullptr, G4ThreeVector(0, 0, zGrease),
fGreaseLV, "Grease_Front", fWrapLV, false, 0);
zPos = zGrease + greaseThickness / 2; // Update position
// Light guide (if used)
if (fUseLightGuide) {
G4double zLightGuide = zPos + fLightGuideLength / 2;
new G4PVPlacement(nullptr, G4ThreeVector(0, 0, zLightGuide),
fLightGuideLV, "LightGuide", fWrapLV, false, 0);
zPos = zLightGuide + fLightGuideLength / 2;
}
// PMT
G4double zPMT = zPos + fPMTLength / 2;
fPMTFrontPV = new G4PVPlacement(nullptr, G4ThreeVector(0, 0, zPMT),
fPMTLV, "PMT", fWrapLV, false, 0);
// Photocathode (at inner surface of PMT closest to scintillator)
G4double zPhotoLocalFront = -fPMTLength / 2 + photoThickness / 2;
fPhotoFrontPV = new G4PVPlacement(nullptr,
G4ThreeVector(0, 0, zPhotoLocalFront),
fPhotoLV,
"Photocathode",
fPMTLV,
false,
0);
// PMT shield
G4double zShield = zPos + (fPMTLength / 2 + shieldThickness / 2);
new G4PVPlacement(nullptr, G4ThreeVector(0, 0, zShield),
shieldLV, "PMTShield", fWrapLV, false, 0);
}
Add Properties
void DetectorConstruction::SetScintMPT(G4Material* m, const G4String& tag) {
std::vector<G4double> E; BuildEnergyGrid(E);
std::vector<G4double> fast(E.size()), slow(E.size()), rindex(E.size()), absl(E.size());
G4double fastTau, slowTau, fracFast, yield, birks;
if (tag == "EJ200") {
// EJ-200: Fast plastic scintillator
for (size_t i = 0; i < E.size(); ++i) {
const G4double Ev = E[i] / eV;
fast[i] = gaussian(Ev, 3.0, 0.3); // Peak at 413 nm
slow[i] = 0.1 * gaussian(Ev, 3.0, 0.4); // Small slow component
rindex[i] = 1.58;
absl[i] = 250. * cm;
}
fastTau = 2.1 * ns;
slowTau = 10.0 * ns;
fracFast = 0.95; // Mostly fast
yield = 10000. / MeV;
birks = 0.126 * mm / MeV;
}
else if (tag == "EJ212") {
// EJ-212: General purpose plastic (similar to EJ-200)
for (size_t i = 0; i < E.size(); ++i) {
const G4double Ev = E[i] / eV;
fast[i] = gaussian(Ev, 2.95, 0.32); // Slightly different emission
slow[i] = 0.15 * gaussian(Ev, 2.95, 0.45);
rindex[i] = 1.58;
absl[i] = 220. * cm;
}
fastTau = 2.4 * ns;
slowTau = 12.0 * ns;
fracFast = 0.90;
yield = 9500. / MeV;
birks = 0.130 * mm / MeV;
}
else if (tag == "EJ276") {
// EJ-276: PSD plastic
for (size_t i = 0; i < E.size(); ++i) {
const G4double Ev = E[i] / eV;
fast[i] = 0.7 * gaussian(Ev, 2.95, 0.25) + 0.3 * gaussian(Ev, 2.7, 0.2);
slow[i] = 0.5 * gaussian(Ev, 2.85, 0.3) + 0.5 * gaussian(Ev, 2.6, 0.25);
rindex[i] = 1.58;
absl[i] = 200. * cm;
}
fastTau = 3.0 * ns;
slowTau = 17.0 * ns;
fracFast = 0.30; // Strong PSD capability
yield = 9000. / MeV;
birks = 0.115 * mm / MeV;
}
else if (tag == "EJ309") {
// EJ-309: Liquid scintillator
for (size_t i = 0; i < E.size(); ++i) {
const G4double Ev = E[i] / eV;
fast[i] = gaussian(Ev, 2.95, 0.35); // Peak at 420 nm
slow[i] = 0.8 * gaussian(Ev, 2.95, 0.4); // Strong slow component
rindex[i] = 1.505;
absl[i] = 300. * cm; // Very clear liquid
}
fastTau = 3.7 * ns;
slowTau = 32.0 * ns;
fracFast = 0.23;
yield = 12000. / MeV;
birks = 0.148 * mm / MeV;
}
else {
// Default - fallback to EJ-200 properties
for (size_t i = 0; i < E.size(); ++i) {
const G4double Ev = E[i] / eV;
fast[i] = gaussian(Ev, 3.0, 0.3);
slow[i] = 0.2 * gaussian(Ev, 3.0, 0.4);
rindex[i] = 1.58;
absl[i] = 200. * cm;
}
fastTau = 2.0 * ns;
slowTau = 10.0 * ns;
fracFast = 0.8;
yield = 8000. / MeV;
birks = 0.1 * mm / MeV;
}
auto* mpt = new G4MaterialPropertiesTable();
mpt->AddProperty("RINDEX", E.data(), rindex.data(), (G4int)E.size(), true);
mpt->AddProperty("ABSLENGTH", E.data(), absl.data(), (G4int)E.size(), true);
mpt->AddProperty("FASTCOMPONENT", E.data(), fast.data(), (G4int)E.size(), true);
mpt->AddProperty("SLOWCOMPONENT", E.data(), slow.data(), (G4int)E.size(), true);
mpt->AddConstProperty("SCINTILLATIONYIELD", yield, true);
mpt->AddConstProperty("RESOLUTIONSCALE", 1.0, true);
mpt->AddConstProperty("FASTTIMECONSTANT", fastTau, true);
mpt->AddConstProperty("SLOWTIMECONSTANT", slowTau, true);
mpt->AddConstProperty("YIELDRATIO", fracFast, true);
// Legacy aliases for compatibility
mpt->AddProperty("SCINTILLATIONCOMPONENT1", E.data(), fast.data(), (G4int)E.size(), true);
mpt->AddProperty("SCINTILLATIONCOMPONENT2", E.data(), slow.data(), (G4int)E.size(), true);
mpt->AddConstProperty("SCINTILLATIONTIMECONSTANT1", fastTau, true);
mpt->AddConstProperty("SCINTILLATIONTIMECONSTANT2", slowTau, true);
mpt->AddConstProperty("SCINTILLATIONYIELD1", fracFast, true);
mpt->AddConstProperty("SCINTILLATIONYIELD2", 1.0 - fracFast, true);
// Realistic timing
mpt->AddConstProperty("RISETIME1", 0.5 * ns, true);
mpt->AddConstProperty("RISETIME2", 1.0 * ns, true);
// Particle-dependent yields (important for PSD scintillators)
mpt->AddConstProperty("ELECTRONSCINTILLATIONYIELD", yield, true);
mpt->AddConstProperty("POSITRONSCINTILLATIONYIELD", yield, true);
mpt->AddConstProperty("PROTONSCINTILLATIONYIELD", 0.65 * yield, true);
mpt->AddConstProperty("DEUTERONSCINTILLATIONYIELD", 0.70 * yield, true);
mpt->AddConstProperty("TRITONSCINTILLATIONYIELD", 0.72 * yield, true);
mpt->AddConstProperty("ALPHASCINTILLATIONYIELD", 0.50 * yield, true);
mpt->AddConstProperty("IONSCINTILLATIONYIELD", 0.60 * yield, true);
m->SetMaterialPropertiesTable(mpt);
m->GetIonisation()->SetBirksConstant(birks);
}
// --------------- construct geometry ---------------
G4VPhysicalVolume* DetectorConstruction::Construct() {
// Clean previous geometry
G4GeometryManager::GetInstance()->OpenGeometry();
G4PhysicalVolumeStore::GetInstance()->Clean();
G4LogicalVolumeStore::GetInstance()->Clean();
G4SolidStore::GetInstance()->Clean();
// Clear existing visualization attributes
for (auto visAtt : fVisAttributes) {
delete visAtt;
}
fVisAttributes.clear();
// Build geometry
SetupVisAttributes();
SetupGeometry();
SetupOpticalSurfaces();
// Geometry validation
G4cout << G4endl;
G4cout << "=== GEOMETRY VALIDATION ===" << G4endl;
G4cout << "Checking for overlaps (tolerance: 0.1 mm)..." << G4endl;
bool overlapCheck = fWorldPV->CheckOverlaps(1000, 0.1 * mm, false);
if (overlapCheck) {
G4cerr << "WARNING: Overlaps detected in geometry!" << G4endl;
}
// Print geometry summary
G4cout << "=== GEOMETRY SUMMARY ===" << G4endl;
G4cout << "Scintillator: " << fScSize / mm << " mm cube of " << fScintMat->GetName() << G4endl;
G4cout << "PMT Type: " << fPMTType << " (" << 2 * fPMTRadius / mm << " mm diameter)" << G4endl;
G4cout << "Wrapping: " << fWrapType << " (" << fWrapReflectivity * 100 << "% reflectivity)" << G4endl;
G4cout << "Light Guide: " << (fUseLightGuide ? "YES" : "NO");
if (fUseLightGuide) G4cout << " (" << fLightGuideLength / mm << " mm length)";
G4cout << G4endl;
G4cout << "=========================" << G4endl;
if (fDumpMPTAtInit) DumpAllMPTs(false);
return fWorldPV;
}
Physicslist
PhysicsList::PhysicsList()
: G4VModularPhysicsList(),
fMessenger(nullptr),
fEmPhysicsList(nullptr),
fOptical(nullptr),
fEmName("standard"),
fScintillationEnabled(true),
fCerenkovEnabled(false),
fWLSEnabled(false),
fRayleighEnabled(false),
fAbsorptionEnabled(true),
fBoundaryEnabled(true)
{
SetVerboseLevel(1);
fMessenger = new PhysicsListMessenger(this);
// Default EM physics
fEmName = "standard";
fEmPhysicsList = new PhysListEmStandard(fEmName);
RegisterPhysics(fEmPhysicsList);
// Core physics
RegisterPhysics(new G4EmExtraPhysics());
RegisterPhysics(new G4DecayPhysics());
RegisterPhysics(new G4RadioactiveDecayPhysics());
RegisterPhysics(new G4HadronPhysicsQGSP_BIC_HP());
// Optional neutron tracking cut
auto* ncut = new G4NeutronTrackingCut();
ncut->SetTimeLimit(NEUTRON_TIME_LIMIT);
RegisterPhysics(ncut);
// Optical physics
fOptical = new G4OpticalPhysics();
RegisterPhysics(fOptical);
ConfigureEMParameters();
ConfigureOpticalParameters();
}
// ---------------------------------------------------------
PhysicsList::~PhysicsList()
{
delete fMessenger;
}
// ---------------------------------------------------------
void PhysicsList::ConstructParticle()
{
G4BosonConstructor().ConstructParticle();
G4LeptonConstructor().ConstructParticle();
G4MesonConstructor().ConstructParticle();
G4BaryonConstructor().ConstructParticle();
G4IonConstructor().ConstructParticle();
G4ShortLivedConstructor().ConstructParticle();
}
// ---------------------------------------------------------
void PhysicsList::ConstructProcess()
{
AddTransportation();
G4VModularPhysicsList::ConstructProcess();
}
// ---------------------------------------------------------
void PhysicsList::SetCuts()
{
SetDefaultCutValue(DEFAULT_CUT_VALUE);
G4VUserPhysicsList::SetCuts();
G4ProductionCutsTable::GetProductionCutsTable()->SetEnergyRange(100 * eV, 1 * GeV);
if (verboseLevel > 0) DumpCutValuesTable();
}
// ---------------------------------------------------------
void PhysicsList::AddPhysicsList(const G4String& name)
{
if (name == fEmName) return;
G4VPhysicsConstructor* newList = nullptr;
if (name == "standard")
newList = new PhysListEmStandard(name);
else if (name == "livermore")
newList = new PhysListEmLivermore(name);
else if (name == "penelope")
newList = new PhysListEmPenelope(name);
else {
G4cout << "[PhysicsList] Unknown EM physics: " << name << G4endl;
return;
}
fEmPhysicsList = newList;
fEmName = name;
ReplacePhysics(fEmPhysicsList);
G4RunManager::GetRunManager()->PhysicsHasBeenModified();
G4cout << "[PhysicsList] Switched EM physics to: " << fEmName << G4endl;
}
// ---------------------------------------------------------
void PhysicsList::ConfigureEMParameters()
{
auto* ep = G4EmParameters::Instance();
ep->SetDefaults();
ep->SetVerbose(verboseLevel);
ep->SetMinEnergy(100 * eV);
ep->SetMaxEnergy(10 * TeV);
ep->SetNumberOfBinsPerDecade(20);
ep->SetMscRangeFactor(0.02);
ep->SetFluo(true);
ep->SetAuger(true);
ep->SetPixe(true);
ep->SetBuildCSDARange(true);
ep->SetLossFluctuations(true);
ep->SetIntegral(false);
}
// ---------------------------------------------------------
void PhysicsList::ConfigureOpticalParameters()
{
auto* op = G4OpticalParameters::Instance();
op->SetProcessActivation("Scintillation", fScintillationEnabled);
op->SetProcessActivation("Cerenkov", fCerenkovEnabled);
op->SetProcessActivation("OpAbsorption", fAbsorptionEnabled);
op->SetProcessActivation("OpBoundary", fBoundaryEnabled);
op->SetProcessActivation("OpRayleigh", fRayleighEnabled);
op->SetProcessActivation("WLS", fWLSEnabled);
op->SetScintTrackSecondariesFirst(true);
op->SetScintStackPhotons(true);
op->SetScintByParticleType(false);
}
// ---------------------------------------------------------
void PhysicsList::SetProcessActive(const G4String& pname, G4bool on)
{
auto* op = G4OpticalParameters::Instance();
op->SetProcessActivation(pname, on);
G4RunManager::GetRunManager()->PhysicsHasBeenModified();
}
// ---------------------------------------------------------
G4bool PhysicsList::IsProcessActive(const G4String& pname) const
{
auto* op = G4OpticalParameters::Instance();
if (pname == "Scintillation") return fScintillationEnabled;
if (pname == "Cerenkov") return fCerenkovEnabled;
if (pname == "OpAbsorption") return fAbsorptionEnabled;
if (pname == "OpBoundary") return fBoundaryEnabled;
if (pname == "OpRayleigh") return fRayleighEnabled;
if (pname == "WLS") return fWLSEnabled;
return false;
}
// ---------------------------------------------------------
void PhysicsList::ListActiveProcesses() const
{
G4cout << "=== Active Optical Processes ===" << G4endl;
G4cout << "Scintillation: " << (fScintillationEnabled ? "ON" : "OFF") << G4endl;
G4cout << "Cerenkov: " << (fCerenkovEnabled ? "ON" : "OFF") << G4endl;
G4cout << "OpAbsorption: " << (fAbsorptionEnabled ? "ON" : "OFF") << G4endl;
G4cout << "OpBoundary: " << (fBoundaryEnabled ? "ON" : "OFF") << G4endl;
G4cout << "OpRayleigh: " << (fRayleighEnabled ? "ON" : "OFF") << G4endl;
G4cout << "WLS: " << (fWLSEnabled ? "ON" : "OFF") << G4endl;
G4cout << "==============================" << G4endl;
}