Skip to content

Commit

Permalink
fix shadowing parameter for distance calculation and move relocation …
Browse files Browse the repository at this point in the history 
…such that NextState is correct when calling RecordHit (#2)

Co-authored-by: Severin Diederichs <severin.diederichs@cern.ch>
  • Loading branch information
@SeverinDiederichs @agheata
2 people authored and agheata committed Sep 13, 2024
1 parent bffd743 commit edb9c39
Show file tree
Hide file tree
Showing 3 changed files with 70 additions and 56 deletions.
36 changes: 22 additions & 14 deletions include/AdePT/kernels/electrons.cuh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -170,20 +170,21 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
// also need to carry them over!

// Check if there's a volume boundary in between.
bool propagated = true;
bool propagated = true;
long hitsurf_index = -1;
double geometryStepLength;
vecgeom::NavigationState nextState;
if (BzFieldValue != 0) {
geometryStepLength = fieldPropagatorBz.ComputeStepAndNextVolume<AdePTNavigator>(
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;
}
Expand Down Expand Up @@ -250,6 +251,18 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
localTime += deltaTime;
properTime += deltaTime * (restMass / eKin);

if (nextState.IsOnBoundary()) {
// if the particle hit a boundary, and is neither stopped or outside, relocate to have the correct next state
// before RecordHit is called
if (!stopped && !nextState.IsOutside()) {
#ifdef ADEPT_USE_SURF
AdePTNavigator::RelocateToNextVolume(pos, dir, hitsurf_index, nextState);
#else
AdePTNavigator::RelocateToNextVolume(pos, dir, nextState);
#endif
}
}

if (auxData.fSensIndex >= 0)
adept_scoring::RecordHit(userScoring, currentTrack.parentID,
IsElectron ? 0 : 1, // Particle type
Expand Down Expand Up @@ -312,19 +325,14 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr

// Kill the particle if it left the world.
if (!nextState.IsOutside()) {
#ifdef ADEPT_USE_SURF
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;
// 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();
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)
Expand Down
20 changes: 11 additions & 9 deletions include/AdePT/kernels/gammas.cuh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -95,16 +95,18 @@ __global__ void TransportGammas(adept::TrackManager<Track> *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
Expand All @@ -130,15 +132,15 @@ __global__ void TransportGammas(adept::TrackManager<Track> *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
Expand Down
70 changes: 37 additions & 33 deletions src/AdePTGeant4Integration.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -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);
}
}
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -423,25 +424,24 @@ 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
aPreStepPoint->SetMaterialCutsCouple(aPreG4TouchableHandle->GetVolume()->GetLogicalVolume()->GetMaterialCutsCouple());
// 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
Expand All @@ -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
Expand Down

0 comments on commit edb9c39

Please sign in to comment.