Skip to content

Commit 55df05f

Browse files
AlvaroEzqpre-commit-ci[bot]
AlvaroEzq
and
pre-commit-ci[bot]
authored
Differentiate between identical volumes (#142)
* filter sensitive physical volumes * save the step info with the correct particular physical volume * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * store physical volumes rotation in world frame * fix primary generator for nested volumes and fix rotation in GetPositionOnGDMLVolume and GetPositionOnGDMLSurface * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * fix world volume as sensitive and active volume * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * exclude daughters volumes from primary generation * save also the rotation into TRestGeant4PrimaryGeneratorInfo * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
1 parent e5ef0ff commit 55df05f

File tree

5 files changed

+141
-22
lines changed

5 files changed

+141
-22
lines changed

‎include/DetectorConstruction.h‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -17,14 +17,19 @@ class DetectorConstruction : public G4VUserDetectorConstruction {
1717

1818
G4GDMLParser* fGdmlParser;
1919
G4VSolid* fGeneratorSolid;
20+
G4LogicalVolume* fGeneratorLogicalVolume;
2021
G4ThreeVector fGeneratorTranslation;
22+
G4RotationMatrix fGeneratorRotation;
2123

2224
Double_t fBoundBoxXMin, fBoundBoxXMax, fBoundBoxYMin, fBoundBoxYMax, fBoundBoxZMin, fBoundBoxZMax;
2325

2426
public:
2527
G4VPhysicalVolume* GetPhysicalVolume(const G4String& physVolName) const;
28+
bool IsPointInsideAnyDaughterVolume(const G4LogicalVolume* logVol, const G4ThreeVector& point) const;
2629
inline G4VSolid* GetGeneratorSolid() const { return fGeneratorSolid; }
30+
inline G4LogicalVolume* GetGeneratorLogicalVolume() const { return fGeneratorLogicalVolume; }
2731
inline G4ThreeVector GetGeneratorTranslation() const { return fGeneratorTranslation; }
32+
inline G4RotationMatrix GetGeneratorRotation() const { return fGeneratorRotation; }
2833

2934
inline Double_t GetBoundBoxXMin() const { return fBoundBoxXMin; }
3035
inline Double_t GetBoundBoxXMax() const { return fBoundBoxXMax; }

‎src/DataModel.cxx‎

Lines changed: 19 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -178,8 +178,25 @@ void TRestGeant4Hits::InsertStep(const G4Step* step) {
178178

179179
const auto& geometryInfo = metadata->GetGeant4GeometryInfo();
180180

181-
const auto& volumeNameGeant4 = step->GetPreStepPoint()->GetPhysicalVolume()->GetName();
182-
const auto& volumeName = geometryInfo.GetAlternativeNameFromGeant4PhysicalName(volumeNameGeant4);
181+
// Get the full name (path) of the physical volume which uniquely identifies it
182+
auto th = step->GetPreStepPoint()->GetTouchable();
183+
G4int depth = th->GetHistoryDepth();
184+
G4String geant4path = "";
185+
if (depth == 0) { // it is the world volume
186+
geant4path = th->GetVolume()->GetName();
187+
}
188+
for (G4int i = 1; i <= depth; ++i) { // start from 1 to skip world volume
189+
// Move the touchable to level i (0 = current volume, depth = world)
190+
G4VPhysicalVolume* pv = th->GetVolume(depth - i);
191+
if (pv) {
192+
if (geant4path != "") {
193+
geant4path += geometryInfo.GetPathSeparator().Data();
194+
}
195+
geant4path += pv->GetName();
196+
}
197+
}
198+
// convert to the names used in gdml (due to assemblies)
199+
const auto volumeName = geometryInfo.GetAlternativePathFromGeant4Path(geant4path);
183200

184201
if (!metadata->IsActiveVolume(volumeName) && step->GetTrack()->GetCurrentStepNumber() != 0) {
185202
// we always store the first step

‎src/DetectorConstruction.cxx‎

Lines changed: 67 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -169,14 +169,26 @@ G4VPhysicalVolume* DetectorConstruction::Construct() {
169169
exit(1);
170170
}
171171

172-
fGeneratorTranslation = gdmlPhysicalVolume->GetTranslation();
172+
TVector3 genTranslation = geometryInfo.GetPosition(
173+
primaryGeneratorInfo.GetSpatialGeneratorFrom().Data()); // in world coordinates
174+
fGeneratorTranslation = {genTranslation.x(), genTranslation.y(), genTranslation.z()};
175+
TRotation genRotation = geometryInfo.GetRotation(
176+
primaryGeneratorInfo.GetSpatialGeneratorFrom().Data()); // in world coordinates
177+
double angle;
178+
TVector3 axis;
179+
genRotation.AngleAxis(angle, axis);
180+
fGeneratorRotation = G4RotationMatrix(G4ThreeVector(axis.X(), axis.Y(), axis.Z()), angle);
173181
if (spatialGeneratorTypeEnum == TRestGeant4PrimaryGeneratorTypes::SpatialGeneratorTypes::SURFACE ||
174182
spatialGeneratorTypeEnum == TRestGeant4PrimaryGeneratorTypes::SpatialGeneratorTypes::VOLUME) {
175183
restG4Metadata->fGeant4PrimaryGeneratorInfo.fSpatialGeneratorPosition = {
176184
fGeneratorTranslation.x(), fGeneratorTranslation.y(), fGeneratorTranslation.z()};
185+
restG4Metadata->fGeant4PrimaryGeneratorInfo.fSpatialGeneratorRotationAxis = {axis.X(), axis.Y(),
186+
axis.Z()};
187+
restG4Metadata->fGeant4PrimaryGeneratorInfo.fSpatialGeneratorRotationValue = angle;
177188
}
178189

179190
fGeneratorSolid = gdmlPhysicalVolume->GetLogicalVolume()->GetSolid();
191+
fGeneratorLogicalVolume = gdmlPhysicalVolume->GetLogicalVolume();
180192

181193
fBoundBoxXMax = -1.e30;
182194
fBoundBoxYMax = -1.e30;
@@ -314,26 +326,64 @@ void DetectorConstruction::ConstructSDandField() {
314326
}
315327
}
316328

329+
bool DetectorConstruction::IsPointInsideAnyDaughterVolume(const G4LogicalVolume* logVol,
330+
const G4ThreeVector& point) const {
331+
if (!logVol) {
332+
G4cout << "DetectorConstruction::IsPointInsideAnyDaughterVolume : logVol is nullptr" << G4endl;
333+
return false;
334+
}
335+
336+
for (size_t i = 0; i < logVol->GetNoDaughters(); ++i) {
337+
G4VPhysicalVolume* daughter = logVol->GetDaughter(i);
338+
G4LogicalVolume* daughterLogVol = daughter->GetLogicalVolume();
339+
G4VSolid* daughterSolid = daughterLogVol->GetSolid();
340+
341+
// localPoint will be the point in the daughter volume reference system
342+
// which is the one used by the solid to determine if the point is inside or not
343+
G4ThreeVector localPoint = point - daughter->GetTranslation();
344+
if (daughter->GetRotation()) {
345+
localPoint = daughter->GetRotation()->inverse() * localPoint;
346+
}
347+
348+
if (daughterSolid->Inside(localPoint) == kInside) {
349+
return true;
350+
}
351+
}
352+
return false;
353+
}
354+
317355
void TRestGeant4GeometryInfo::PopulateFromGeant4World(const G4VPhysicalVolume* world) {
318356
auto detector = (DetectorConstruction*)G4RunManager::GetRunManager()->GetUserDetectorConstruction();
319357
TRestGeant4Metadata* restG4Metadata = detector->fSimulationManager->GetRestMetadata();
320358

321359
// Recursive function to traverse the nested volume geometry
322-
std::function<void(const G4VPhysicalVolume*, size_t&, const G4String pathSoFar, const G4ThreeVector&)>
360+
std::function<void(const G4VPhysicalVolume*, size_t&, const G4String pathSoFar, const G4ThreeVector&,
361+
const G4RotationMatrix&)>
323362
ProcessVolumeRecursively = [&](const G4VPhysicalVolume* volume, size_t& index,
324-
const G4String pathSoFar, const G4ThreeVector& parentPosition) {
363+
const G4String pathSoFar, const G4ThreeVector& parentPosition,
364+
const G4RotationMatrix& parentRotation) {
325365
G4String currentPath = pathSoFar;
326366
if (volume->GetName() != world->GetName()) { // avoid all paths including 'world_PV/' at the
327367
// beginning
328368
currentPath += (currentPath.empty() ? "" : fPathSeparator.Data()) + volume->GetName();
329369
}
330-
G4ThreeVector positionInWorld = volume->GetTranslation() + parentPosition;
370+
G4ThreeVector localPosition = volume->GetTranslation();
371+
localPosition = parentRotation * localPosition;
372+
G4RotationMatrix localRotation =
373+
volume->GetRotation() ? *volume->GetRotation() : G4RotationMatrix(); // identity if nullptr
374+
375+
// accumulate the position
376+
G4ThreeVector positionInWorld = parentPosition + localPosition;
377+
378+
// accumulate the rotation
379+
G4RotationMatrix rotationInWorld = parentRotation;
380+
rotationInWorld *= localRotation;
331381

332382
// First process the daughters to have the same order in volume IDs as before
333383
G4LogicalVolume* logVol = volume->GetLogicalVolume();
334384
for (size_t i = 0; i < logVol->GetNoDaughters(); ++i) {
335385
G4VPhysicalVolume* daughter = logVol->GetDaughter(i);
336-
ProcessVolumeRecursively(daughter, index, currentPath, positionInWorld);
386+
ProcessVolumeRecursively(daughter, index, currentPath, positionInWorld, rotationInWorld);
337387
}
338388

339389
// Process this volume
@@ -351,6 +401,14 @@ void TRestGeant4GeometryInfo::PopulateFromGeant4World(const G4VPhysicalVolume* w
351401
fLogicalToPhysicalMap[nameLogical].emplace_back(namePhysical);
352402
fPhysicalToPositionInWorldMap[physicalNewName] = {positionInWorld.x(), positionInWorld.y(),
353403
positionInWorld.z()};
404+
// Convert G4RotationMatrix to TRotation and store it
405+
double angle;
406+
G4ThreeVector axis;
407+
rotationInWorld.getAngleAxis(angle, axis);
408+
TRotation rotInWorld;
409+
rotInWorld.Rotate(angle, TVector3(axis.x(), axis.y(), axis.z()));
410+
fPhysicalToRotationInWorldMap[physicalNewName] = rotInWorld;
411+
354412
InsertVolumeName(index, physicalNewName);
355413
/*
356414
std::cout << "Index: " << index << std::endl;
@@ -360,6 +418,9 @@ void TRestGeant4GeometryInfo::PopulateFromGeant4World(const G4VPhysicalVolume* w
360418
std::cout << "\tMaterial: " << nameMaterial << std::endl;
361419
std::cout << "\tPosition: (" << positionInWorld.x() << ", " << positionInWorld.y() << ", " <<
362420
positionInWorld.z() << ")mm" << std::endl;
421+
std::cout << "\tRotation (angle, axis): (" << angle << ", (" << axis.x() << ", " << axis.y() << ",
422+
"
423+
<< axis.z() << "))" << std::endl;
363424
*/
364425
if (!fIsAssembly &&
365426
GetAlternativeNameFromGeant4PhysicalName(namePhysical).Data() != namePhysical) {
@@ -378,5 +439,5 @@ void TRestGeant4GeometryInfo::PopulateFromGeant4World(const G4VPhysicalVolume* w
378439

379440
// Start recursion from world volume
380441
size_t index = 0;
381-
ProcessVolumeRecursively(world, index, "", G4ThreeVector(0, 0, 0));
442+
ProcessVolumeRecursively(world, index, "", G4ThreeVector(0, 0, 0), G4RotationMatrix());
382443
}

‎src/PrimaryGeneratorAction.cxx‎

Lines changed: 15 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -751,11 +751,17 @@ void PrimaryGeneratorAction::GenPositionOnGDMLVolume(double& x, double& y, doubl
751751
x = xMin + (xMax - xMin) * G4UniformRand();
752752
y = yMin + (yMax - yMin) * G4UniformRand();
753753
z = zMin + (zMax - zMin) * G4UniformRand();
754-
} while (detector->GetGeneratorSolid()->Inside(G4ThreeVector(x, y, z)) != kInside);
754+
} while (detector->GetGeneratorSolid()->Inside(G4ThreeVector(x, y, z)) != kInside ||
755+
detector->IsPointInsideAnyDaughterVolume(detector->GetGeneratorLogicalVolume(),
756+
G4ThreeVector(x, y, z)));
755757

756-
x = x + detector->GetGeneratorTranslation().x();
757-
y = y + detector->GetGeneratorTranslation().y();
758-
z = z + detector->GetGeneratorTranslation().z();
758+
// Rotate and translate since the solid has no sense of the physical volume rotation or translation
759+
G4ThreeVector position = G4ThreeVector(x, y, z);
760+
position = detector->GetGeneratorRotation() * position;
761+
762+
x = position.x() + detector->GetGeneratorTranslation().x();
763+
y = position.y() + detector->GetGeneratorTranslation().y();
764+
z = position.z() + detector->GetGeneratorTranslation().z();
759765
}
760766

761767
void PrimaryGeneratorAction::GenPositionOnGDMLSurface(double& x, double& y, double& z) {
@@ -766,13 +772,12 @@ void PrimaryGeneratorAction::GenPositionOnGDMLSurface(double& x, double& y, doub
766772

767773
G4ThreeVector position = detector->GetGeneratorSolid()->GetPointOnSurface();
768774

769-
x = position.x();
770-
y = position.y();
771-
z = position.z();
775+
// Rotate and translate since the solid has no sense of the physical volume rotation or translation
776+
position = detector->GetGeneratorRotation() * position;
772777

773-
x = x + detector->GetGeneratorTranslation().x();
774-
y = y + detector->GetGeneratorTranslation().y();
775-
z = z + detector->GetGeneratorTranslation().z();
778+
x = position.x() + detector->GetGeneratorTranslation().x();
779+
y = position.y() + detector->GetGeneratorTranslation().y();
780+
z = position.z() + detector->GetGeneratorTranslation().z();
776781
}
777782

778783
void PrimaryGeneratorAction::GenPositionOnBoxVolume(double& x, double& y, double& z) {

‎src/SensitiveDetector.cxx‎

Lines changed: 35 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -16,10 +16,41 @@ SensitiveDetector::SensitiveDetector(SimulationManager* simulationManager, const
1616

1717
G4bool SensitiveDetector::ProcessHits(G4Step* step, G4TouchableHistory*) {
1818
// return value will always be ignored, its present for backwards compatibility (I guess)
19-
const auto volumeName = fSimulationManager->GetRestMetadata()
20-
->GetGeant4GeometryInfo()
21-
.GetAlternativeNameFromGeant4PhysicalName(
22-
step->GetPreStepPoint()->GetPhysicalVolume()->GetName());
19+
TRestGeant4Metadata* restG4Metadata = fSimulationManager->GetRestMetadata();
20+
const auto& geometryInfo = restG4Metadata->GetGeant4GeometryInfo();
21+
const auto volumeName = geometryInfo.GetAlternativeNameFromGeant4PhysicalName(
22+
step->GetPreStepPoint()->GetPhysicalVolume()->GetName());
23+
24+
/* Sensitive detector are assigned to logical volumes, but several physical volumes
25+
can share the same logical volume. Now we get the physical volume where the pre-step
26+
is located and filter if it is one of the declared sensitive volumes by the user*/
27+
28+
// Get the full name (path) of the physical volume which uniquely identifies it
29+
auto th = step->GetPreStepPoint()->GetTouchable();
30+
G4int depth = th->GetHistoryDepth();
31+
G4String geant4path = "";
32+
if (depth == 0) { // it is the world volume
33+
geant4path = th->GetVolume()->GetName();
34+
}
35+
for (G4int i = 1; i <= depth; ++i) { // start from 1 to skip world volume
36+
// Move the touchable to level i (0 = current volume, depth = world)
37+
G4VPhysicalVolume* pv = th->GetVolume(depth - i);
38+
if (pv) {
39+
if (geant4path != "") {
40+
geant4path += geometryInfo.GetPathSeparator().Data();
41+
}
42+
geant4path += pv->GetName();
43+
}
44+
}
45+
// convert to the names used in gdml (due to assemblies)
46+
const auto physicalFullName = geometryInfo.GetAlternativePathFromGeant4Path(geant4path);
47+
48+
// search and filter if this physical volume is sensitive
49+
auto sensitiveVolumes = restG4Metadata->GetSensitiveVolumes();
50+
if (std::find(sensitiveVolumes.begin(), sensitiveVolumes.end(), physicalFullName) ==
51+
sensitiveVolumes.end()) {
52+
return false;
53+
}
2354

2455
const bool isGeantino = step->GetTrack()->GetParticleDefinition() == G4Geantino::Definition();
2556

0 commit comments

Comments
 (0)