Skip to content

Detector z positions #13

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
tobias-boeckh merged 4 commits into from
Dec 9, 2025
Merged

Detector z positions #13

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
b1f1eba
Consider scintillator thickness in calculation of start z position
tobias-boeckh Dec 9, 2025
217d46f
Fix calculation of z vertex position
tobias-boeckh Dec 9, 2025
8301692
Add flag to set silicon box thickness
tobias-boeckh Dec 9, 2025
df74fa3
Add silicon box thickness to readme and fix unit
tobias-boeckh Dec 9, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 13 additions & 9 deletions Pinpoint/include/DetectorConstruction.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,15 +33,16 @@ 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 SetBoxThickness(G4double thickness) { fBoxThickness = thickness; }
void SetNLayers(G4int nLayers) { fNLayers = nLayers; }
void SetPixelHeight(G4double height) { fPixelHeight = height; }
void SetPixelWidth(G4double width) { fPixelWidth = width; }
Expand Down Expand Up @@ -75,21 +76,22 @@ class DetectorConstruction : public G4VUserDetectorConstruction
std::vector<G4double> GetPixelZPositions() const
{
std::vector<G4double> 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;
};

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; }
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) {
Expand All @@ -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;
Expand All @@ -122,6 +125,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;

Expand Down
1 change: 1 addition & 0 deletions Pinpoint/include/DetectorConstructionMessenger.hh
View file
Open in desktop
Copy link
Contributor

@benw22022 benw22022 Dec 9, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

New geometry commands should be documented in README table

Copy link
Contributor Author

@tobias-boeckh tobias-boeckh Dec 9, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @benw22022, I've updated the readme and fixed the unit.

Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@ class DetectorConstructionMessenger: public G4UImessenger {

G4UIcmdWithADoubleAndUnit* tungstenThicknessCmd;
G4UIcmdWithADoubleAndUnit* siliconThicknessCmd;
G4UIcmdWithADoubleAndUnit* boxThicknessCmd;
G4UIcmdWithAnInteger* nLayersCmd;
G4UIcmdWithADoubleAndUnit* pixelHeightCmd;
G4UIcmdWithADoubleAndUnit* pixelWidthCmd;
Expand Down
3 changes: 3 additions & 0 deletions Pinpoint/macros/test.mac
View file
Open in desktop
Copy link
Contributor

@benw22022 benw22022 Dec 9, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

should say # default category is mm

Original file line number Diff line number Diff line change
Expand Up @@ -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 mm

# Number of sampling layers
/det/setNLayers 100

Expand Down
31 changes: 15 additions & 16 deletions Pinpoint/src/DetectorConstruction.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -174,18 +173,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<int>(100.0*cm / layerThickness); // maximum layers allowed
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<int>(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;
Expand Down Expand Up @@ -221,18 +220,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());
Expand All @@ -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
Expand Down Expand Up @@ -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 * 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 * 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);
Expand Down
11 changes: 11 additions & 0 deletions Pinpoint/src/DetectorConstructionMessenger.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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");
Expand Down Expand Up @@ -104,6 +110,7 @@ DetectorConstructionMessenger::~DetectorConstructionMessenger() {
delete detDir;
delete tungstenThicknessCmd;
delete siliconThicknessCmd;
delete boxThicknessCmd;
delete nLayersCmd;
delete pixelHeightCmd;
delete pixelWidthCmd;
Expand All @@ -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);
Expand Down
25 changes: 5 additions & 20 deletions Pinpoint/src/generators/GFaserGenerator.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
}

1 change: 1 addition & 0 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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` |
Expand Down