From edb9c391861604e0d6d37c2785f416ebf33ef3f0 Mon Sep 17 00:00:00 2001 From: Severin Diederichs <65728274+SeverinDiederichs@users.noreply.github.com> Date: Thu, 12 Sep 2024 11:33:39 +0200 Subject: [PATCH] fix shadowing parameter for distance calculation and move relocation such that NextState is correct when calling RecordHit (#2) Co-authored-by: Severin Diederichs --- include/AdePT/kernels/electrons.cuh | 36 +++++++++------ include/AdePT/kernels/gammas.cuh | 20 +++++---- src/AdePTGeant4Integration.cpp | 70 +++++++++++++++-------------- 3 files changed, 70 insertions(+), 56 deletions(-) diff --git a/include/AdePT/kernels/electrons.cuh b/include/AdePT/kernels/electrons.cuh index 82c93e31..4440c541 100644 --- a/include/AdePT/kernels/electrons.cuh +++ b/include/AdePT/kernels/electrons.cuh @@ -170,20 +170,21 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager( - eKin, restMass, Charge, geometricalStepLengthFromPhysics, pos, dir, navState, nextState, hitsurf_index, propagated, safety); + eKin, restMass, Charge, geometricalStepLengthFromPhysics, pos, dir, navState, nextState, hitsurf_index, + propagated, safety); } else { #ifdef ADEPT_USE_SURF - double geometryStepLength = AdePTNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics, - navState, nextState, hitsurf_index, kPush); + geometryStepLength = AdePTNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics, + navState, nextState, hitsurf_index, kPush); #else - double geometryStepLength = AdePTNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics, - navState, nextState, kPush); + geometryStepLength = AdePTNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics, + navState, nextState, kPush); #endif pos += geometryStepLength * dir; } @@ -250,6 +251,18 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager= 0) adept_scoring::RecordHit(userScoring, currentTrack.parentID, IsElectron ? 0 : 1, // Particle type @@ -312,19 +325,14 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManagerGetLogicalVolume()->id(); + const int nextlvolID = navState.Top()->GetLogicalVolume()->id(); #else - const int nextlvolID = navState.GetLogicalId(); + const int nextlvolID = navState.GetLogicalId(); #endif VolAuxData const &nextauxData = auxDataArray[nextlvolID]; if (nextauxData.fGPUregion > 0) diff --git a/include/AdePT/kernels/gammas.cuh b/include/AdePT/kernels/gammas.cuh index e6fcc2db..d5db4b31 100644 --- a/include/AdePT/kernels/gammas.cuh +++ b/include/AdePT/kernels/gammas.cuh @@ -95,16 +95,18 @@ __global__ void TransportGammas(adept::TrackManager *gammas, Secondaries // Check if there's a volume boundary in between. vecgeom::NavigationState nextState; + double geometryStepLength; #ifdef ADEPT_USE_SURF long hitsurf_index = -1; - double geometryStepLength = AdePTNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics, - navState, nextState, hitsurf_index, kPush); + geometryStepLength = AdePTNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics, navState, + nextState, hitsurf_index, kPush); #else - double geometryStepLength = AdePTNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics, - navState, nextState, kPush); + geometryStepLength = AdePTNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics, navState, + nextState, kPush); #endif - // printf("pvol=%d step=%g onboundary=%d pos={%g, %g, %g} dir={%g, %g, %g}\n", navState.TopId(), geometryStepLength, - // nextState.IsOnBoundary(), pos[0], pos[1], pos[2], dir[0], dir[1], dir[2]); + // printf("pvol=%d step=%g onboundary=%d pos={%g, %g, %g} dir={%g, %g, %g}\n", navState.TopId(), + // geometryStepLength, + // nextState.IsOnBoundary(), pos[0], pos[1], pos[2], dir[0], dir[1], dir[2]); pos += geometryStepLength * dir; // Set boundary state in navState so the next step and secondaries get the @@ -130,15 +132,15 @@ __global__ void TransportGammas(adept::TrackManager *gammas, Secondaries // Kill the particle if it left the world. if (!nextState.IsOutside()) { #ifdef ADEPT_USE_SURF - AdePTNavigator::RelocateToNextVolume(pos, dir, hitsurf_index, nextState); + AdePTNavigator::RelocateToNextVolume(pos, dir, hitsurf_index, nextState); if (nextState.IsOutside()) continue; #else AdePTNavigator::RelocateToNextVolume(pos, dir, nextState); #endif // Move to the next boundary. navState = nextState; - //printf(" -> pvol=%d pos={%g, %g, %g} \n", navState.TopId(), pos[0], pos[1], pos[2]); - // Check if the next volume belongs to the GPU region and push it to the appropriate queue + // printf(" -> pvol=%d pos={%g, %g, %g} \n", navState.TopId(), pos[0], pos[1], pos[2]); + // Check if the next volume belongs to the GPU region and push it to the appropriate queue #ifndef ADEPT_USE_SURF const int nextlvolID = navState.Top()->GetLogicalVolume()->id(); #else diff --git a/src/AdePTGeant4Integration.cpp b/src/AdePTGeant4Integration.cpp index 1133efc4..372c7330 100644 --- a/src/AdePTGeant4Integration.cpp +++ b/src/AdePTGeant4Integration.cpp @@ -278,17 +278,20 @@ void AdePTGeant4Integration::ProcessGPUHits(HostScoring &aScoring, HostScoring:: // Reconstruct G4NavigationHistory and G4Step, and call the SD code for each hit for (size_t i = aStats.fBufferStart; i < aStats.fBufferStart + aStats.fUsedSlots; i++) { // Get Hit index (Circular buffer) - int aHitIdx = i % aScoring.fBufferCapacity; - - vecgeom::NavigationState &aNavState = aScoring.fGPUHitsBuffer_host[aHitIdx].fPreStepPoint.fNavigationState; + int aHitIdx = i % aScoring.fBufferCapacity; + vecgeom::NavigationState &preNavState = aScoring.fGPUHitsBuffer_host[aHitIdx].fPreStepPoint.fNavigationState; // Reconstruct Pre-Step point G4NavigationHistory - FillG4NavigationHistory(aNavState, fPreG4NavigationHistory); + FillG4NavigationHistory(preNavState, fPreG4NavigationHistory); ((G4TouchableHistory *)fPreG4TouchableHistoryHandle()) ->UpdateYourself(fPreG4NavigationHistory->GetTopVolume(), fPreG4NavigationHistory); // Reconstruct Post-Step point G4NavigationHistory - FillG4NavigationHistory(aNavState, fPostG4NavigationHistory); - ((G4TouchableHistory *)fPostG4TouchableHistoryHandle()) - ->UpdateYourself(fPostG4NavigationHistory->GetTopVolume(), fPostG4NavigationHistory); + vecgeom::NavigationState &postNavState = aScoring.fGPUHitsBuffer_host[aHitIdx].fPostStepPoint.fNavigationState; + if (!postNavState + .IsOutside()) { // if it is outside, don't set anything and check for nullptr GetVolume in FillG4Step + FillG4NavigationHistory(postNavState, fPostG4NavigationHistory); + ((G4TouchableHistory *)fPostG4TouchableHistoryHandle()) + ->UpdateYourself(fPostG4NavigationHistory->GetTopVolume(), fPostG4NavigationHistory); + } // Reconstruct G4Step switch (aScoring.fGPUHitsBuffer_host[aHitIdx].fParticleType) { @@ -317,7 +320,8 @@ void AdePTGeant4Integration::ProcessGPUHits(HostScoring &aScoring, HostScoring:: } } -void AdePTGeant4Integration::FillG4NavigationHistory(vecgeom::NavigationState aNavState, G4NavigationHistory *aG4NavigationHistory) +void AdePTGeant4Integration::FillG4NavigationHistory(vecgeom::NavigationState aNavState, + G4NavigationHistory *aG4NavigationHistory) { // Get the current depth of the history (corresponding to the previous reconstructed touchable) auto aG4HistoryDepth = aG4NavigationHistory->GetDepth(); @@ -351,13 +355,10 @@ void AdePTGeant4Integration::FillG4NavigationHistory(vecgeom::NavigationState aN aG4HistoryDepth = aLevel; } else { // If the navigation state is deeper than the current history we need to add the new levels - if(aLevel) - { + if (aLevel) { aG4NavigationHistory->NewLevel(pnewvol, kNormal, pnewvol->GetCopyNo()); aG4HistoryDepth++; - } - else - { + } else { aG4NavigationHistory->SetFirstEntry(pnewvol); } } @@ -369,15 +370,15 @@ void AdePTGeant4Integration::FillG4NavigationHistory(vecgeom::NavigationState aN void AdePTGeant4Integration::FillG4Step(GPUHit *aGPUHit, G4Step *aG4Step, G4TouchableHandle &aPreG4TouchableHandle, G4TouchableHandle &aPostG4TouchableHandle) { - const G4ThreeVector *aPostStepPointMomentumDirection = new G4ThreeVector(aGPUHit->fPostStepPoint.fMomentumDirection.x(), - aGPUHit->fPostStepPoint.fMomentumDirection.y(), - aGPUHit->fPostStepPoint.fMomentumDirection.z()); - const G4ThreeVector *aPostStepPointPolarization = new G4ThreeVector(aGPUHit->fPostStepPoint.fPolarization.x(), - aGPUHit->fPostStepPoint.fPolarization.y(), - aGPUHit->fPostStepPoint.fPolarization.z()); - const G4ThreeVector *aPostStepPointPosition = new G4ThreeVector(aGPUHit->fPostStepPoint.fPosition.x(), - aGPUHit->fPostStepPoint.fPosition.y(), - aGPUHit->fPostStepPoint.fPosition.z()); + const G4ThreeVector *aPostStepPointMomentumDirection = + new G4ThreeVector(aGPUHit->fPostStepPoint.fMomentumDirection.x(), aGPUHit->fPostStepPoint.fMomentumDirection.y(), + aGPUHit->fPostStepPoint.fMomentumDirection.z()); + const G4ThreeVector *aPostStepPointPolarization = + new G4ThreeVector(aGPUHit->fPostStepPoint.fPolarization.x(), aGPUHit->fPostStepPoint.fPolarization.y(), + aGPUHit->fPostStepPoint.fPolarization.z()); + const G4ThreeVector *aPostStepPointPosition = + new G4ThreeVector(aGPUHit->fPostStepPoint.fPosition.x(), aGPUHit->fPostStepPoint.fPosition.y(), + aGPUHit->fPostStepPoint.fPosition.z()); // G4Step aG4Step->SetStepLength(aGPUHit->fStepLength); // Real data @@ -391,8 +392,8 @@ void AdePTGeant4Integration::FillG4Step(GPUHit *aGPUHit, G4Step *aG4Step, G4Touc // G4Track G4Track *aTrack = aG4Step->GetTrack(); - aTrack->SetTrackID(aGPUHit->fParentID); // Missing data - aTrack->SetParentID(aGPUHit->fParentID); // ID of the initial particle that entered AdePT + aTrack->SetTrackID(aGPUHit->fParentID); // Missing data + aTrack->SetParentID(aGPUHit->fParentID); // ID of the initial particle that entered AdePT aTrack->SetPosition(*aPostStepPointPosition); // Real data // aTrack->SetGlobalTime(0); // Missing data // aTrack->SetLocalTime(0); // Missing data @@ -423,16 +424,15 @@ void AdePTGeant4Integration::FillG4Step(GPUHit *aGPUHit, G4Step *aG4Step, G4Touc // Pre-Step Point G4StepPoint *aPreStepPoint = aG4Step->GetPreStepPoint(); - aPreStepPoint->SetPosition(G4ThreeVector(aGPUHit->fPreStepPoint.fPosition.x(), - aGPUHit->fPreStepPoint.fPosition.y(), - aGPUHit->fPreStepPoint.fPosition.z())); // Real data + aPreStepPoint->SetPosition(G4ThreeVector(aGPUHit->fPreStepPoint.fPosition.x(), aGPUHit->fPreStepPoint.fPosition.y(), + aGPUHit->fPreStepPoint.fPosition.z())); // Real data // aPreStepPoint->SetLocalTime(0); // Missing data // aPreStepPoint->SetGlobalTime(0); // Missing data // aPreStepPoint->SetProperTime(0); // Missing data aPreStepPoint->SetMomentumDirection(G4ThreeVector(aGPUHit->fPreStepPoint.fMomentumDirection.x(), aGPUHit->fPreStepPoint.fMomentumDirection.y(), aGPUHit->fPreStepPoint.fMomentumDirection.z())); // Real data - aPreStepPoint->SetKineticEnergy(aGPUHit->fPreStepPoint.fEKin); // Real data + aPreStepPoint->SetKineticEnergy(aGPUHit->fPreStepPoint.fEKin); // Real data // aPreStepPoint->SetVelocity(0); // Missing data aPreStepPoint->SetTouchableHandle(aPreG4TouchableHandle); // Real data aPreStepPoint->SetMaterial(aPreG4TouchableHandle->GetVolume()->GetLogicalVolume()->GetMaterial()); // Real data @@ -440,8 +440,8 @@ void AdePTGeant4Integration::FillG4Step(GPUHit *aGPUHit, G4Step *aG4Step, G4Touc // aPreStepPoint->SetSensitiveDetector(nullptr); // Missing data // aPreStepPoint->SetSafety(0); // Missing data aPreStepPoint->SetPolarization(G4ThreeVector(aGPUHit->fPreStepPoint.fPolarization.x(), - aGPUHit->fPreStepPoint.fPolarization.y(), - aGPUHit->fPreStepPoint.fPolarization.z())); // Real data + aGPUHit->fPreStepPoint.fPolarization.y(), + aGPUHit->fPreStepPoint.fPolarization.z())); // Real data // aPreStepPoint->SetStepStatus(G4StepStatus::fUndefined); // Missing data // aPreStepPoint->SetProcessDefinedStep(nullptr); // Missing data // aPreStepPoint->SetMass(0); // Missing data @@ -458,9 +458,13 @@ void AdePTGeant4Integration::FillG4Step(GPUHit *aGPUHit, G4Step *aG4Step, G4Touc aPostStepPoint->SetMomentumDirection(*aPostStepPointMomentumDirection); // Real data aPostStepPoint->SetKineticEnergy(aGPUHit->fPostStepPoint.fEKin); // Real data // aPostStepPoint->SetVelocity(0); // Missing data - aPostStepPoint->SetTouchableHandle(aPostG4TouchableHandle); // Real data - aPostStepPoint->SetMaterial(aPostG4TouchableHandle->GetVolume()->GetLogicalVolume()->GetMaterial()); // Real data - aPostStepPoint->SetMaterialCutsCouple(aPostG4TouchableHandle->GetVolume()->GetLogicalVolume()->GetMaterialCutsCouple()); + if (fPostG4TouchableHistoryHandle->GetVolume()) { // protect against nullptr if postNavState is outside + aPostStepPoint->SetTouchableHandle(fPostG4TouchableHistoryHandle); // Real data + aPostStepPoint->SetMaterial( + fPostG4TouchableHistoryHandle->GetVolume()->GetLogicalVolume()->GetMaterial()); // Real data + aPostStepPoint->SetMaterialCutsCouple( + fPostG4TouchableHistoryHandle->GetVolume()->GetLogicalVolume()->GetMaterialCutsCouple()); + } // aPostStepPoint->SetSensitiveDetector(nullptr); // Missing data // aPostStepPoint->SetSafety(0); // Missing data aPostStepPoint->SetPolarization(*aPostStepPointPolarization); // Real data