From b1f1eba42879461ad1976e92b3eab0ce02ce374e Mon Sep 17 00:00:00 2001 From: Tobias Boeckh Date: Tue, 9 Dec 2025 11:11:02 +0100 Subject: [PATCH 1/4] Consider scintillator thickness in calculation of start z position --- Pinpoint/include/DetectorConstruction.hh | 19 ++++++++++--------- Pinpoint/src/DetectorConstruction.cc | 24 ++++++++++++------------ 2 files changed, 22 insertions(+), 21 deletions(-) diff --git a/Pinpoint/include/DetectorConstruction.hh b/Pinpoint/include/DetectorConstruction.hh index f0dcacd..fbc3b22 100644 --- a/Pinpoint/include/DetectorConstruction.hh +++ b/Pinpoint/include/DetectorConstruction.hh @@ -33,14 +33,14 @@ class DetectorConstruction : public G4VUserDetectorConstruction if (thickness <= 0) { G4cerr << "Error: Tungsten thickness must be positive." << G4endl; return; - } + } - fTungstenThickness = thickness; - G4cout << "Set tungsten thickness to " << fTungstenThickness/mm << " mm" << G4endl; + fTungstenThickness = thickness; + G4cout << "Set tungsten thickness to " << fTungstenThickness/mm << " mm" << G4endl; - // Optional: trigger geometry rebuild - // G4RunManager::GetRunManager()->ReinitializeGeometry(); -} + // Optional: trigger geometry rebuild + // G4RunManager::GetRunManager()->ReinitializeGeometry(); + } void SetSiliconThickness(G4double thickness) { fSiliconThickness = thickness; } void SetNLayers(G4int nLayers) { fNLayers = nLayers; } void SetPixelHeight(G4double height) { fPixelHeight = height; } @@ -75,10 +75,9 @@ class DetectorConstruction : public G4VUserDetectorConstruction std::vector GetPixelZPositions() const { std::vector zPositions; - G4double layerSpacing = fTungstenThickness + fSiliconThickness; - G4double startZ = -((fNLayers - 1) * layerSpacing) / 2; + G4double startZ = -0.5 * fNLayers * fLayerThickness; for (G4int i = 0; i < fNLayers; ++i) { - zPositions.push_back(startZ + i * layerSpacing); + zPositions.push_back(startZ + i * fLayerThickness + fTungstenThickness + fSiliconThickness / 2); } return zPositions; }; @@ -90,6 +89,7 @@ class DetectorConstruction : public G4VUserDetectorConstruction G4double GetPixelWidth() const { return fPixelWidth; } G4double GetDetectorWidth() const { return fDetectorWidth; } G4double GetDetectorHeight() const { return fDetectorHeight; } + G4double GetLayerThickness() const { return fLayerThickness; } G4int GetSimFlag() const {return sim_flag;} G4int GetScintBarFlag() const { if (false) { @@ -122,6 +122,7 @@ class DetectorConstruction : public G4VUserDetectorConstruction G4double fScintBarWidth = 9.85 * mm; G4double fScintBarHeight = 9.80 * mm; G4double fScintThickness = 5.0 * mm; + G4double fLayerThickness = 0.0 * mm; G4int sim_flag = 0; G4bool scint_bar_flag = true; diff --git a/Pinpoint/src/DetectorConstruction.cc b/Pinpoint/src/DetectorConstruction.cc index 128dbb1..02adbc9 100644 --- a/Pinpoint/src/DetectorConstruction.cc +++ b/Pinpoint/src/DetectorConstruction.cc @@ -174,18 +174,18 @@ G4VPhysicalVolume* DetectorConstruction::Construct() G4cout << "Silicon thickness per layer: " << fSiliconThickness/um << " um" << G4endl; G4cout << "Creating " << nPixelsX << " x " << nPixelsY << " pixels per silicon layer" << G4endl; G4cout << "Pixel size: " << fPixelWidth/micrometer << " x " << fPixelHeight/micrometer << " μm" << G4endl; - G4double layerThickness = fTungstenThickness + boxThickness + fSiliconThickness; - if(sim_flag == -1) { layerThickness = fTungstenThickness + boxThickness + fSiliconThickness;} //only pixel, TPTPTPTP... - if(sim_flag == 0) { layerThickness = 2.0*fTungstenThickness + boxThickness + fSiliconThickness + fScintThickness;} //pixel + single scintillator, TPTSTPTS... - if(sim_flag == 1) { layerThickness = 2.0*fTungstenThickness + boxThickness + fSiliconThickness + 2.0*fScintThickness;} //pixel + double scintillator, TPTSSTPTSS... - //auto layerThickness = fTungstenThickness + boxThickness + fSiliconThickness; - auto maxLayers = static_cast(100.0*cm / layerThickness); // maximum layers allowed + fLayerThickness = fTungstenThickness + boxThickness + fSiliconThickness; + if(sim_flag == -1) { fLayerThickness = fTungstenThickness + boxThickness + fSiliconThickness;} //only pixel, TPTPTPTP... + if(sim_flag == 0) { fLayerThickness = 2.0*fTungstenThickness + boxThickness + fSiliconThickness + fScintThickness;} //pixel + single scintillator, TPTSTPTS... + if(sim_flag == 1) { fLayerThickness = 2.0*fTungstenThickness + boxThickness + fSiliconThickness + 2.0*fScintThickness;} //pixel + double scintillator, TPTSSTPTSS... + //auto fLayerThickness = fTungstenThickness + boxThickness + fSiliconThickness; + auto maxLayers = static_cast(100.0*cm / fLayerThickness); // maximum layers allowed if(fNLayers > maxLayers) { G4cout << "Warning: Reducing number of layers from " << fNLayers << " to " << maxLayers << " to keep detector thickness <= 100 cm." << G4endl; fNLayers = maxLayers; } - auto detectorThickness = fNLayers * layerThickness; + auto detectorThickness = fNLayers * fLayerThickness; auto worldSizeX = 1.2 * fDetectorWidth; auto worldSizeY = 1.2 * fDetectorHeight; auto worldSizeZ = 1.2 * detectorThickness; @@ -221,18 +221,18 @@ G4VPhysicalVolume* DetectorConstruction::Construct() new G4PVPlacement(0, G4ThreeVector(0., 0., 0.5 * detectorThickness), detectorLV, "Detector", worldLV, false, 0, fCheckOverlaps); // Layer - auto layerS = new G4Box("Layer", 0.5 * fDetectorWidth, 0.5 * fDetectorHeight, layerThickness / 2); + auto layerS = new G4Box("Layer", 0.5 * fDetectorWidth, 0.5 * fDetectorHeight, fLayerThickness / 2); auto layerLV = new G4LogicalVolume(layerS, worldMaterial, "Layer"); - fLayerPV = new G4PVReplica("Layer", layerLV, detectorLV, kZAxis, fNLayers, layerThickness); + fLayerPV = new G4PVReplica("Layer", layerLV, detectorLV, kZAxis, fNLayers, fLayerThickness); //Cursor for placements inside each layer - G4double zCursor = -0.5*layerThickness; + G4double zCursor = -0.5*fLayerThickness; // Tungsten auto tungstenS = new G4Box("Tungsten", 0.5 * fDetectorWidth, 0.5 * fDetectorHeight, 0.5 * fTungstenThickness); auto tungstenLV = new G4LogicalVolume(tungstenS, tungstenMaterial, "Tungsten"); zCursor += 0.5*fTungstenThickness; - //fTarget_phys.push_back(new G4PVPlacement(0, G4ThreeVector(0., 0., -0.5 * layerThickness + 0.5 * fTungstenThickness), tungstenLV, "Tungsten", layerLV, false, 0, fCheckOverlaps)); + //fTarget_phys.push_back(new G4PVPlacement(0, G4ThreeVector(0., 0., -0.5 * fLayerThickness + 0.5 * fTungstenThickness), tungstenLV, "Tungsten", layerLV, false, 0, fCheckOverlaps)); new G4PVPlacement(0, G4ThreeVector(0., 0., zCursor), tungstenLV, "Tungsten", layerLV, false, 0, fCheckOverlaps); G4VisAttributes* TargetVisAtt = new G4VisAttributes(G4Colour::Red()); @@ -368,7 +368,7 @@ if(sim_flag == 1) // // Box // auto boxS = new G4Box("Box", 0.5 * fDetectorWidth, 0.5 * fDetectorHeight, 0.5 * boxThickness); // auto boxLV = new G4LogicalVolume(boxS, worldMaterial, "Box"); - // new G4PVPlacement(nullptr, G4ThreeVector(0., 0., -0.5 * layerThickness + fTungstenThickness + fSiliconThickness + 0.5 * boxThickness), boxLV, "Box", layerLV, false, 0, fCheckOverlaps); + // new G4PVPlacement(nullptr, G4ThreeVector(0., 0., -0.5 * fLayerThickness + fTungstenThickness + fSiliconThickness + 0.5 * boxThickness), boxLV, "Box", layerLV, false, 0, fCheckOverlaps); // G4cout << "Detector consists of " << fNLayers << " layers of: [ " << fTungstenThickness / mm << "mm of " << tungstenMaterial->GetName() << " + " << fSiliconThickness / mm << "mm of " // << siliconMaterial->GetName() << " + " << boxThickness / mm << "mm of " << worldMaterial->GetName() << " ] " << G4endl; From 217d46f96e3f71f9ba7a78f821b81b13eb38ad11 Mon Sep 17 00:00:00 2001 From: Tobias Boeckh Date: Tue, 9 Dec 2025 11:11:52 +0100 Subject: [PATCH 2/4] Fix calculation of z vertex position --- Pinpoint/src/generators/GFaserGenerator.cc | 25 +++++----------------- 1 file changed, 5 insertions(+), 20 deletions(-) diff --git a/Pinpoint/src/generators/GFaserGenerator.cc b/Pinpoint/src/generators/GFaserGenerator.cc index e6aefc0..b6ac719 100644 --- a/Pinpoint/src/generators/GFaserGenerator.cc +++ b/Pinpoint/src/generators/GFaserGenerator.cc @@ -231,27 +231,12 @@ G4bool GFaserGenerator::FindParticleDefinition(G4int const pdg, G4ParticleDefini G4double GFaserGenerator::GenerateRandomZVertex(G4int layerIndex) const { auto *runManager = G4RunManager::GetRunManager(); auto detector = (DetectorConstruction*) (runManager->GetUserDetectorConstruction()); - G4VPhysicalVolume* layerPV = detector->GetLayerPhysVol(); - if (!layerPV) { - G4cerr << "Error: Layer volume not found!" << G4endl; - return 0.0; - } - - EAxis axis; - G4int nReplicas; - G4double width; - G4double offset; - G4bool consuming; - layerPV->GetReplicationData(axis, nReplicas, width, offset, consuming); - - if (layerIndex >= nReplicas) { - G4cerr << "Error: Layer index " << layerIndex << " is out of range." - << "Only " << nReplicas << " layers available." << G4endl; - return 0.0; - } - - G4double z = (layerIndex + G4UniformRand()) * width + offset; + G4double layerThickness = detector->GetLayerThickness(); + G4double tungstenThickness = detector->GetTungstenThickness(); + G4int nLayers = detector->GetNumberOfLayers(); + G4double startZ = -0.5 * nLayers * layerThickness; + G4double z = startZ + layerIndex * layerThickness + tungstenThickness * G4UniformRand(); return z/m; // convert to meters } From 8301692774b665d388a0aee532de35c0768cd85d Mon Sep 17 00:00:00 2001 From: Tobias Boeckh Date: Tue, 9 Dec 2025 11:21:48 +0100 Subject: [PATCH 3/4] Add flag to set silicon box thickness --- Pinpoint/include/DetectorConstruction.hh | 3 +++ .../include/DetectorConstructionMessenger.hh | 1 + Pinpoint/macros/test.mac | 3 +++ Pinpoint/src/DetectorConstruction.cc | 19 +++++++++---------- Pinpoint/src/DetectorConstructionMessenger.cc | 11 +++++++++++ 5 files changed, 27 insertions(+), 10 deletions(-) diff --git a/Pinpoint/include/DetectorConstruction.hh b/Pinpoint/include/DetectorConstruction.hh index fbc3b22..d5adc59 100644 --- a/Pinpoint/include/DetectorConstruction.hh +++ b/Pinpoint/include/DetectorConstruction.hh @@ -42,6 +42,7 @@ class DetectorConstruction : public G4VUserDetectorConstruction // G4RunManager::GetRunManager()->ReinitializeGeometry(); } void SetSiliconThickness(G4double thickness) { fSiliconThickness = thickness; } + void SetBoxThickness(G4double thickness) { fBoxThickness = thickness; } void SetNLayers(G4int nLayers) { fNLayers = nLayers; } void SetPixelHeight(G4double height) { fPixelHeight = height; } void SetPixelWidth(G4double width) { fPixelWidth = width; } @@ -84,6 +85,7 @@ class DetectorConstruction : public G4VUserDetectorConstruction G4double GetTungstenThickness() const { return fTungstenThickness; } G4double GetSiliconThickness() const { return fSiliconThickness; } + G4double GetBoxThickness() const { return fBoxThickness; } G4int GetNumberOfLayers() const { return fNLayers; } G4double GetPixelHeight() const { return fPixelHeight; } G4double GetPixelWidth() const { return fPixelWidth; } @@ -112,6 +114,7 @@ class DetectorConstruction : public G4VUserDetectorConstruction G4double fTungstenThickness = 5 * mm; G4double fSiliconThickness = 50 * um; + G4double fBoxThickness = 5.0 * mm; G4int fNLayers = 100; G4double fPixelHeight = 22.8 * um; G4double fPixelWidth = 20.8 * um; diff --git a/Pinpoint/include/DetectorConstructionMessenger.hh b/Pinpoint/include/DetectorConstructionMessenger.hh index 2da5bee..2c12c5e 100644 --- a/Pinpoint/include/DetectorConstructionMessenger.hh +++ b/Pinpoint/include/DetectorConstructionMessenger.hh @@ -32,6 +32,7 @@ class DetectorConstructionMessenger: public G4UImessenger { G4UIcmdWithADoubleAndUnit* tungstenThicknessCmd; G4UIcmdWithADoubleAndUnit* siliconThicknessCmd; + G4UIcmdWithADoubleAndUnit* boxThicknessCmd; G4UIcmdWithAnInteger* nLayersCmd; G4UIcmdWithADoubleAndUnit* pixelHeightCmd; G4UIcmdWithADoubleAndUnit* pixelWidthCmd; diff --git a/Pinpoint/macros/test.mac b/Pinpoint/macros/test.mac index 3cdca89..55c1ac2 100644 --- a/Pinpoint/macros/test.mac +++ b/Pinpoint/macros/test.mac @@ -14,6 +14,9 @@ # Set silicon sensor thickness /det/setSiliconThickness 50 um # default category is um +# Set silicon box thickness +/det/setBoxThickness 5.0 mm # default category is um + # Number of sampling layers /det/setNLayers 100 diff --git a/Pinpoint/src/DetectorConstruction.cc b/Pinpoint/src/DetectorConstruction.cc index 02adbc9..bf9d2dd 100644 --- a/Pinpoint/src/DetectorConstruction.cc +++ b/Pinpoint/src/DetectorConstruction.cc @@ -158,7 +158,6 @@ G4VPhysicalVolume* DetectorConstruction::Construct() { // Geometry parameters // https://iopscience.iop.org/article/10.1088/1748-0221/20/02/C02015 - G4double boxThickness = fTungstenThickness - fSiliconThickness; G4bool fCheckOverlaps = true; //cleaning scintillator logical volume container @@ -174,11 +173,11 @@ G4VPhysicalVolume* DetectorConstruction::Construct() G4cout << "Silicon thickness per layer: " << fSiliconThickness/um << " um" << G4endl; G4cout << "Creating " << nPixelsX << " x " << nPixelsY << " pixels per silicon layer" << G4endl; G4cout << "Pixel size: " << fPixelWidth/micrometer << " x " << fPixelHeight/micrometer << " μm" << G4endl; - fLayerThickness = fTungstenThickness + boxThickness + fSiliconThickness; - if(sim_flag == -1) { fLayerThickness = fTungstenThickness + boxThickness + fSiliconThickness;} //only pixel, TPTPTPTP... - if(sim_flag == 0) { fLayerThickness = 2.0*fTungstenThickness + boxThickness + fSiliconThickness + fScintThickness;} //pixel + single scintillator, TPTSTPTS... - if(sim_flag == 1) { fLayerThickness = 2.0*fTungstenThickness + boxThickness + fSiliconThickness + 2.0*fScintThickness;} //pixel + double scintillator, TPTSSTPTSS... - //auto fLayerThickness = fTungstenThickness + boxThickness + fSiliconThickness; + fLayerThickness = fTungstenThickness + fBoxThickness + fSiliconThickness; + if(sim_flag == -1) { fLayerThickness = fTungstenThickness + fBoxThickness + fSiliconThickness;} //only pixel, TPTPTPTP... + if(sim_flag == 0) { fLayerThickness = 2.0*fTungstenThickness + fBoxThickness + fSiliconThickness + fScintThickness;} //pixel + single scintillator, TPTSTPTS... + if(sim_flag == 1) { fLayerThickness = 2.0*fTungstenThickness + fBoxThickness + fSiliconThickness + 2.0*fScintThickness;} //pixel + double scintillator, TPTSSTPTSS... + //auto fLayerThickness = fTungstenThickness + fBoxThickness + fSiliconThickness; auto maxLayers = static_cast(100.0*cm / fLayerThickness); // maximum layers allowed if(fNLayers > maxLayers) { G4cout << "Warning: Reducing number of layers from " << fNLayers @@ -252,7 +251,7 @@ G4VPhysicalVolume* DetectorConstruction::Construct() ScintLayerAtrrib->SetVisibility(true); ScintLayerAtrrib->SetForceSolid(true); - zCursor += 0.5*fTungstenThickness + boxThickness + 0.5*fSiliconThickness; + zCursor += 0.5*fTungstenThickness + fBoxThickness + 0.5*fSiliconThickness; auto silicofNLayers = new G4Box("SiliconLayer", 0.5 * fDetectorWidth, 0.5 * fDetectorHeight, 0.5 * fSiliconThickness); auto siliconLayerLV = new G4LogicalVolume(silicofNLayers, siliconMaterial, "SiliconLayer"); // Changed to siliconMaterial @@ -366,12 +365,12 @@ if(sim_flag == 1) } // // Box - // auto boxS = new G4Box("Box", 0.5 * fDetectorWidth, 0.5 * fDetectorHeight, 0.5 * boxThickness); + // auto boxS = new G4Box("Box", 0.5 * fDetectorWidth, 0.5 * fDetectorHeight, 0.5 * fBoxThickness); // auto boxLV = new G4LogicalVolume(boxS, worldMaterial, "Box"); - // new G4PVPlacement(nullptr, G4ThreeVector(0., 0., -0.5 * fLayerThickness + fTungstenThickness + fSiliconThickness + 0.5 * boxThickness), boxLV, "Box", layerLV, false, 0, fCheckOverlaps); + // new G4PVPlacement(nullptr, G4ThreeVector(0., 0., -0.5 * fLayerThickness + fTungstenThickness + fSiliconThickness + 0.5 * fBoxThickness), boxLV, "Box", layerLV, false, 0, fCheckOverlaps); // G4cout << "Detector consists of " << fNLayers << " layers of: [ " << fTungstenThickness / mm << "mm of " << tungstenMaterial->GetName() << " + " << fSiliconThickness / mm << "mm of " - // << siliconMaterial->GetName() << " + " << boxThickness / mm << "mm of " << worldMaterial->GetName() << " ] " << G4endl; + // << siliconMaterial->GetName() << " + " << fBoxThickness / mm << "mm of " << worldMaterial->GetName() << " ] " << G4endl; // Write GDML file if it doesn't exist std::ifstream file(fWriteFile); diff --git a/Pinpoint/src/DetectorConstructionMessenger.cc b/Pinpoint/src/DetectorConstructionMessenger.cc index 96e3be9..42bb9c9 100644 --- a/Pinpoint/src/DetectorConstructionMessenger.cc +++ b/Pinpoint/src/DetectorConstructionMessenger.cc @@ -40,6 +40,12 @@ DetectorConstructionMessenger::DetectorConstructionMessenger(DetectorConstructio siliconThicknessCmd->SetParameterName("SiliconThickness", false); siliconThicknessCmd->SetRange("SiliconThickness>0."); + boxThicknessCmd = new G4UIcmdWithADoubleAndUnit("/det/setBoxThickness", this); + boxThicknessCmd->SetUnitCategory("Length"); + boxThicknessCmd->SetDefaultUnit("mm"); + boxThicknessCmd->SetParameterName("BoxThickness", false); + boxThicknessCmd->SetRange("BoxThickness>0."); + nLayersCmd = new G4UIcmdWithAnInteger("/det/setNLayers", this); nLayersCmd->SetParameterName("NLayers", false); nLayersCmd->SetRange("NLayers>0"); @@ -104,6 +110,7 @@ DetectorConstructionMessenger::~DetectorConstructionMessenger() { delete detDir; delete tungstenThicknessCmd; delete siliconThicknessCmd; + delete boxThicknessCmd; delete nLayersCmd; delete pixelHeightCmd; delete pixelWidthCmd; @@ -126,6 +133,10 @@ void DetectorConstructionMessenger::SetNewValue(G4UIcommand* command, G4String n G4double thickness = siliconThicknessCmd->ConvertToDimensionedDouble(newValues); det->SetSiliconThickness(thickness); } + if (command == boxThicknessCmd) { + G4double thickness = boxThicknessCmd->ConvertToDimensionedDouble(newValues); + det->SetBoxThickness(thickness); + } if (command == nLayersCmd) { G4int nLayers = nLayersCmd->GetNewIntValue(newValues); det->SetNLayers(nLayers); From df74fa37e7268b1e0af5ac5ad1f7193a706bb342 Mon Sep 17 00:00:00 2001 From: Tobias Boeckh Date: Tue, 9 Dec 2025 12:46:52 +0100 Subject: [PATCH 4/4] Add silicon box thickness to readme and fix unit --- Pinpoint/macros/test.mac | 2 +- README.md | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/Pinpoint/macros/test.mac b/Pinpoint/macros/test.mac index 55c1ac2..dc19070 100644 --- a/Pinpoint/macros/test.mac +++ b/Pinpoint/macros/test.mac @@ -15,7 +15,7 @@ /det/setSiliconThickness 50 um # default category is um # Set silicon box thickness -/det/setBoxThickness 5.0 mm # default category is um +/det/setBoxThickness 5.0 mm # default category is mm # Number of sampling layers /det/setNLayers 100 diff --git a/README.md b/README.md index 1c8e45a..e76d2dc 100644 --- a/README.md +++ b/README.md @@ -47,6 +47,7 @@ There are a number of user defined macro commands which can be used to control t |:--|:--|:--| |`/det/setTungstenThickness` | Tungsten sheet thickness in mm | `5 mm`| |`/det/setSiliconThickness` | Silicon layer thickness in $\mu$m | `50 um`| +|`/det/setBoxThickness` | Silicon box (filled with air) thickness in mm | `5 mm`| |`/det/setNLayers` | Number of tungsten/pixel layer | `100`| |`/det/setPixelHeight`| Set the height of an individual pixel in $\mu$m | `22.8 um` | |`/det/setPixelWidth`| Set the width of an individual pixel in $\mu$m | `20.8 um` |