Skip to content

Commit 3ec3d73

Browse files
ggalgocziplexoos
andauthored
fix(gpuraytrace): add multi-component scintillation genstep collection (#272)
GPU genstep collection was only using SCINTILLATIONTIMECONSTANT1 (fast component, τ=7ns). LAr has a 25% slow component (SCINTILLATIONTIMECONSTANT2, τ=1400ns) that was never collected. This caused GPU to produce zero hits with arrival time >100ns while Geant4 showed 20% of hits in the 100-11,000ns range. Scintillation fix (src/GPURaytrace.h): - Read up to 3 scintillation components (SCINTILLATIONTIMECONSTANT1/2/3, SCINTILLATIONYIELD1/2/3) from material properties - Split photon count proportionally across components based on yield fractions - Create separate gensteps for each component with the correct time constant - Last component receives rounding remainder to preserve total photon count exactly Shall be merged after rest of the wavelength shifting PRs --------- Co-authored-by: Dmitri Smirnov <dmixsmi@gmail.com>
1 parent 77ba39c commit 3ec3d73

2 files changed

Lines changed: 46 additions & 3 deletions

File tree

src/GPURaytrace.h

Lines changed: 44 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -531,10 +531,51 @@ struct SteppingAction : G4UserSteppingAction
531531
<< G4endl;
532532
return;
533533
}
534-
G4double SCINTILLATIONTIMECONSTANT1 = MPT->GetConstProperty(kSCINTILLATIONTIMECONSTANT1);
534+
// G4 11.x supports up to 3 scintillation components
535+
const G4int tcKeys[3] = {kSCINTILLATIONTIMECONSTANT1, kSCINTILLATIONTIMECONSTANT2,
536+
kSCINTILLATIONTIMECONSTANT3};
537+
const G4int yieldKeys[3] = {kSCINTILLATIONYIELD1, kSCINTILLATIONYIELD2, kSCINTILLATIONYIELD3};
535538

536-
U4::CollectGenstep_DsG4Scintillation_r4695(aTrack, aStep, fNumPhotons, 1,
537-
SCINTILLATIONTIMECONSTANT1);
539+
G4double tc[3] = {0, 0, 0};
540+
G4double yield[3] = {0, 0, 0};
541+
G4double yieldSum = 0;
542+
G4int nComp = 0;
543+
544+
for (G4int c = 0; c < 3; c++)
545+
{
546+
if (MPT->ConstPropertyExists(tcKeys[c]))
547+
{
548+
tc[c] = MPT->GetConstProperty(tcKeys[c]);
549+
yield[c] = MPT->ConstPropertyExists(yieldKeys[c]) ? MPT->GetConstProperty(yieldKeys[c])
550+
: (c == 0 ? 1.0 : 0.0);
551+
yieldSum += yield[c];
552+
nComp = c + 1;
553+
}
554+
}
555+
556+
if (yieldSum <= 0.0)
557+
{
558+
G4cout << "WARNING: scintillation yields sum to <= 0 for material '"
559+
<< aMaterial->GetName() << "'; falling back to single component." << G4endl;
560+
yield[0] = 1.0;
561+
yieldSum = 1.0;
562+
nComp = 1;
563+
}
564+
565+
G4AutoLock lock(&genstep_mutex);
566+
G4int nRemaining = fNumPhotons;
567+
for (G4int c = 0; c < nComp; c++)
568+
{
569+
G4int nPhotComp;
570+
if (c == nComp - 1)
571+
nPhotComp = nRemaining; // last component gets remainder
572+
else
573+
nPhotComp = static_cast<G4int>(fNumPhotons * yield[c] / yieldSum);
574+
nRemaining -= nPhotComp;
575+
576+
if (nPhotComp > 0)
577+
U4::CollectGenstep_DsG4Scintillation_r4695(aTrack, aStep, nPhotComp, c + 1, tc[c]);
578+
}
538579
}
539580
}
540581
}

tests/geom/opticks_raindrop_with_scintillation.gdml

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,7 @@
1919
<matrix coldim="1" name="FAST_TIME_CONST" values="6.0"/>
2020
<matrix coldim="1" name="SLOW_TIME_CONST" values="35.0"/>
2121
<matrix coldim="1" name="YIELD_RATIO" values="0.8"/>
22+
<matrix coldim="1" name="YIELD_RATIO_SLOW" values="0.2"/>
2223
<matrix coldim="2" name="REEMISSION_PROB" values="1.0*eV 0.0 4.0*eV 0.0"/>
2324
</define>
2425

@@ -58,6 +59,7 @@
5859
<property name="SCINTILLATIONTIMECONSTANT1" ref="FAST_TIME_CONST"/>
5960
<property name="SCINTILLATIONTIMECONSTANT2" ref="SLOW_TIME_CONST"/>
6061
<property name="SCINTILLATIONYIELD1" ref="YIELD_RATIO"/>
62+
<property name="SCINTILLATIONYIELD2" ref="YIELD_RATIO_SLOW"/>
6163
<!-- Old Geant4 10.x / Opticks property names -->
6264
<property name="FASTCOMPONENT" ref="SCINT_SPECTRUM"/>
6365
<property name="SLOWCOMPONENT" ref="SCINT_SPECTRUM"/>

0 commit comments

Comments
 (0)