Generic Biasing with a parallel World for multiple types of particles

weller · June 14, 2021, 11:15pm

Hi everyone,
I am trying to implement generic biasing based on geometry in a parallel world. For this I have used example GB06 and a little GB03 as inspiration, however my particles seem to ignore some of the boundaries…
My strategy was to implement nested boxes in the parallel world, i.e. box in box in world, with increasing importance for increasing depth.

What happens is the following:

Gamma gets created in world volume
Entering box 1: splitting into n (e.g. n=32) identical gammas
Entering box 2: each gamma is splitting into n identical gammas, again, so n*n in total
Scattering on some object → lots of random directions

Now things don’t work out:

Leaving box 2 (thus getting back into box 1) → gammas shoud be killed with 1/n ratio, but 100% survive
Leaving box 1 → getting killed with 1/n ratio. Since leaving box 2 does not reduce the number of gammas, there is on the order of n gammas leaving box 1, although it should only be 1 (plus/minus some from the scattering in step 4).

Particles “see” the boundary between boxes, as there is a Step Point drawn in the GUI. However, most of the tracks should end at the boundary, and only 1/n tracks should continue.

It seems to me that particles cloned inside the second box are not placed into this box. In example GB03 I find this comment:

// Check if step is limited by the geometry: as we attached the biasing operator // to the absorber layer, this volume boundary is the one of the absorber. // (check of current step # of track is inelegant, but is to fix a "feature" // that a cloned track can wrongly be seen in the wrong volume, because of numerical // precision issue. In this case it makes a tiny step, which should be disregarded).

Alternatively I implemented G4ImportanceBiasing, but this only works for a single type of particle…

Is there anyone willing to help me work this out? I would then gladly create working code that I can share.
Thank you so much in advance

edit: I think this could be solved by forcing the clones to be in the same volume as the splitted particles? is there a way to do this?

weller · June 16, 2021, 2:12pm

edit: this post is probably not correct, see next post in this thread

my understanding of the problem now after some deeper digging: when the splits/clones are created,

for ( G4int iSplit = 1 ; iSplit <  splittingFactor ; iSplit++ ) {
      G4Track* clone = new G4Track( *track );
      clone->SetWeight( weightOfTrack );
      fParticleChange.AddSecondary( clone );
}

they are placed into the most “senior” mother volume in the parallel world instead of the “youngest” daughter volume. (I double checked what happens if I displace the clones a small amount to be off the geometric boundaries

clone->SetPosition(track->GetPosition()+10.*CLHEP::um*track->GetMomentumDirection());

but no luck )
As clones have already passed the boundary to the youngest daughter volume they will stay in their current volume until crossing another boundary (i.e., leaving the mother volume into the world) where splitting/killing again is applied. the boundary out of the youngest daughter volume however is always skipped…

any ideas how to force the navigator to accept clones into the desired volume? is this a bug?

weller · June 25, 2021, 3:30pm

I have created a working example here:

can’t get it to work so far, but additional insight: it’s not clones getting recognized in the wrong depth of volume…

Now there is an “insideBox” inside an “outsideBox” within the “parallelWorld”, nothing in the mass-world. Particle #1 travels through as expected:

G4WT0 > * G4Track Information:   Particle = gamma,   Track ID = 1,   Parent ID = 0
G4WT0 > *********************************************************************************************************
G4WT0 > 
G4WT0 > Step#    X(mm)    Y(mm)    Z(mm) KinE(MeV)  dE(MeV) StepLeng TrackLeng  NextVolume ProcName
G4WT0 >     0        0        0      200       0.1        0        0         0       World initStep
G4WT0 > 1 / 1 DistanceToApply: parallelWorld -> parallelWorld Current Volume: parallelWorld
G4WT0 > 1 / 1 GenerateBiasingFinalState: parallelWorld -> parallelWorld Current Volume: outsideBox Pos: (0,0,100)
G4WT0 > 1 from id: 2  to id: 3 parallelWorld -> outsideBox
G4WT0 >     1        0        0      100       0.1        0      100       100       World biasLimiter
G4WT0 >     :----- List of 2ndaries - #SpawnInStep=  1(Rest= 0,Along= 0,Post= 1), #SpawnTotal=  1 ---------------
G4WT0 >     :         0         0       100       0.1              gamma
G4WT0 >     :----------------------------------------------------------------- EndOf2ndaries Info ---------------
G4WT0 > 1 / 2 DistanceToApply: parallelWorld -> outsideBox Current Volume: outsideBox
G4WT0 > 1 / 2 GenerateBiasingFinalState: parallelWorld -> outsideBox Current Volume: insideBox Pos: (0,0,50)
G4WT0 > 1 from id: 3  to id: 4 outsideBox -> insideBox
G4WT0 >     2        0        0       50       0.1        0       50       150       World biasLimiter
G4WT0 >     :----- List of 2ndaries - #SpawnInStep=  2(Rest= 0,Along= 0,Post= 2), #SpawnTotal=  3 ---------------
G4WT0 >     :         0         0        50       0.1              gamma
G4WT0 >     :         0         0        50       0.1              gamma
G4WT0 >     :----------------------------------------------------------------- EndOf2ndaries Info ---------------
G4WT0 > 1 / 3 DistanceToApply: outsideBox -> insideBox Current Volume: insideBox
G4WT0 > 1 / 3 GenerateBiasingFinalState: outsideBox -> insideBox Current Volume: outsideBox Pos: (0,0,-50)
G4WT0 > 1 from id: 4  to id: 3 insideBox -> outsideBox
G4WT0 > 1 killed

getting killed at the boundary insideBox/outsideBox.

One of the clones, #2, starts at the boundary parallelWorld/outsideBox, splits when entering insideBox, survives leaving insideBox, gets killed when leaving outsideBox:

G4WT0 > *********************************************************************************************************
G4WT0 > * G4Track Information:   Particle = gamma,   Track ID = 2,   Parent ID = 1
G4WT0 > *********************************************************************************************************
G4WT0 > 
G4WT0 > Step#    X(mm)    Y(mm)    Z(mm) KinE(MeV)  dE(MeV) StepLeng TrackLeng  NextVolume ProcName
G4WT0 >     0        0        0      100       0.1        0        0         0       World initStep
G4WT0 > 2 / 1 DistanceToApply: outsideBox -> outsideBox Current Volume: outsideBox
G4WT0 > 2 / 1 GenerateBiasingFinalState: outsideBox -> outsideBox Current Volume: insideBox Pos: (0,0,50)
G4WT0 > 2 from id: 3  to id: 4 outsideBox -> insideBox
G4WT0 >     1        0        0       50       0.1        0       50        50       World biasLimiter
G4WT0 >     :----- List of 2ndaries - #SpawnInStep=  2(Rest= 0,Along= 0,Post= 2), #SpawnTotal=  2 ---------------
G4WT0 >     :         0         0        50       0.1              gamma
G4WT0 >     :         0         0        50       0.1              gamma
G4WT0 >     :----------------------------------------------------------------- EndOf2ndaries Info ---------------
G4WT0 > 2 / 2 DistanceToApply: outsideBox -> insideBox Current Volume: insideBox
G4WT0 > 2 / 2 GenerateBiasingFinalState: outsideBox -> insideBox Current Volume: outsideBox Pos: (0,0,-50)
G4WT0 > 2 from id: 4  to id: 3 insideBox -> outsideBox
G4WT0 > 2 killed

However, when placing a scattering object into the mass world, tracing particle #1 again:

G4WT0 > *********************************************************************************************************
G4WT0 > * G4Track Information:   Particle = gamma,   Track ID = 1,   Parent ID = 0
G4WT0 > *********************************************************************************************************
G4WT0 > 
G4WT0 > Step#    X(mm)    Y(mm)    Z(mm) KinE(MeV)  dE(MeV) StepLeng TrackLeng  NextVolume ProcName
G4WT0 >     0        0        0      200       0.1        0        0         0       World initStep
G4WT0 > 1 / 1 DistanceToApply: parallelWorld -> parallelWorld Current Volume: parallelWorld
G4WT0 > 1 / 1 GenerateBiasingFinalState: parallelWorld -> parallelWorld Current Volume: outsideBox Pos: (0,0,100)
G4WT0 > 1 from id: 2  to id: 3 parallelWorld -> outsideBox
G4WT0 >     1        0        0      100       0.1        0      100       100       World biasLimiter
G4WT0 >     :----- List of 2ndaries - #SpawnInStep=  1(Rest= 0,Along= 0,Post= 1), #SpawnTotal=  1 ---------------
G4WT0 >     :         0         0       100       0.1              gamma
G4WT0 >     :----------------------------------------------------------------- EndOf2ndaries Info ---------------
G4WT0 > 1 / 2 DistanceToApply: parallelWorld -> outsideBox Current Volume: outsideBox
G4WT0 > 1 / 2 GenerateBiasingFinalState: parallelWorld -> outsideBox Current Volume: insideBox Pos: (0,0,50)
G4WT0 > 1 from id: 3  to id: 4 outsideBox -> insideBox
G4WT0 >     2        0        0       50       0.1        0       50       150       World biasLimiter
G4WT0 >     :----- List of 2ndaries - #SpawnInStep=  2(Rest= 0,Along= 0,Post= 2), #SpawnTotal=  3 ---------------
G4WT0 >     :         0         0        50       0.1              gamma
G4WT0 >     :         0         0        50       0.1              gamma
G4WT0 >     :----------------------------------------------------------------- EndOf2ndaries Info ---------------
G4WT0 > 1 / 3 DistanceToApply: outsideBox -> insideBox Current Volume: insideBox
G4WT0 > 1 / 3 GenerateBiasingFinalState: outsideBox -> insideBox Current Volume: insideBox Pos: (0,0,5)
G4WT0 > 1 from id: 4  to id: 4 insideBox -> insideBox
G4WT0 > 1 not doing anything
G4WT0 >     3        0        0        5       0.1        0       45       195         Box CoupledTransportation
G4WT0 >     4        0        0    -4.58    0.0972 0.000119     9.58       205         Box compt
G4WT0 >     :----- List of 2ndaries - #SpawnInStep=  1(Rest= 0,Along= 0,Post= 1), #SpawnTotal=  4 ---------------
G4WT0 >     :         0         0     -4.58   0.00266                 e-
G4WT0 >     :----------------------------------------------------------------- EndOf2ndaries Info ---------------
G4WT0 >     5   -0.274     0.14       -5    0.0972        0    0.525       205       World CoupledTransportation
G4WT0 >     6    -29.3       15      -50    0.0972        0     55.6       261       World biasLimiter
G4WT0 > 1 / 7 DistanceToApply: insideBox -> outsideBox Current Volume: outsideBox
G4WT0 > 1 / 7 GenerateBiasingFinalState: insideBox -> outsideBox Current Volume: parallelWorld Pos: (-61.6426,31.5805,-100)
G4WT0 > 1 from id: 3  to id: 2 outsideBox -> parallelWorld
G4WT0 >     7    -61.6     31.6     -100    0.0972        0     61.8       322       World biasLimiter
G4WT0 > 1 / 8 DistanceToApply: outsideBox -> parallelWorld Current Volume: parallelWorld
G4WT0 >     8     -159     81.2     -250    0.0972        0      185       508  OutOfWorld CoupledTransportation

There is no transition from insideBox to outsideBox, or id: 4 to id: 3. From the pattern, I would expect this to happen before/at step 6…

I would be so very happy for suggestions on how to debug this

weller · June 29, 2021, 8:39pm

Is it that the biasing is triggered (call of “DistanceToApply” followed by “GenerateBiasingFinalState”) when entering an attached volume?

I have added the positions to the output from above, last block.
Boundary for outer box at z=100mm and inner box at 50mm:

G4WT0 > *********************************************************************************************************
G4WT0 > * G4Track Information:   Particle = gamma,   Track ID = 1,   Parent ID = 0
G4WT0 > *********************************************************************************************************
G4WT0 > 
G4WT0 > Step#    X(mm)    Y(mm)    Z(mm) KinE(MeV)  dE(MeV) StepLeng TrackLeng  NextVolume ProcName
G4WT0 >     0        0        0      200       0.1        0        0         0       World initStep
G4WT0 > 1 / 1 DistanceToApply: parallelWorld -> parallelWorld Current Volume: parallelWorld Pos: (0,0,200)
G4WT0 > 1 / 1 GenerateBiasingFinalState: parallelWorld -> parallelWorld Current Volume: outsideBox Pos: (0,0,200) -> (0,0,100)
G4WT0 > 1 from id: 2  to id: 3 parallelWorld -> outsideBox
G4WT0 > 1 splitting: 2
G4WT0 >     1        0        0      100       0.1        0      100       100       World biasLimiter
G4WT0 >     :----- List of 2ndaries - #SpawnInStep=  1(Rest= 0,Along= 0,Post= 1), #SpawnTotal=  1 ---------------
G4WT0 >     :         0        10       100       0.1              gamma
G4WT0 >     :----------------------------------------------------------------- EndOf2ndaries Info ---------------
G4WT0 > 1 / 2 DistanceToApply: parallelWorld -> outsideBox Current Volume: outsideBox Pos: (0,0,100)
G4WT0 > 1 / 2 GenerateBiasingFinalState: parallelWorld -> outsideBox Current Volume: insideBox Pos: (0,0,100) -> (0,0,50)
G4WT0 > 1 from id: 3  to id: 4 outsideBox -> insideBox
G4WT0 > 1 splitting: 3
G4WT0 >     2        0        0       50       0.1        0       50       150       World biasLimiter
G4WT0 >     :----- List of 2ndaries - #SpawnInStep=  2(Rest= 0,Along= 0,Post= 2), #SpawnTotal=  3 ---------------
G4WT0 >     :         0        10        50       0.1              gamma
G4WT0 >     :        10        10        50       0.1              gamma
G4WT0 >     :----------------------------------------------------------------- EndOf2ndaries Info ---------------
# 
# I think this is the relevant part:
# begin
G4WT0 > 1 / 3 DistanceToApply: outsideBox -> insideBox Current Volume: insideBox Pos: (0,0,50)
G4WT0 > 1 / 3 GenerateBiasingFinalState: outsideBox -> insideBox Current Volume: insideBox Pos: (0,0,50) -> (0,0,5)
G4WT0 > 1 from id: 4  to id: 4 insideBox -> insideBox
G4WT0 > 1 not doing anything
# end
G4WT0 >     3        0        0        5       0.1        0       45       195         Box CoupledTransportation
G4WT0 >     4        0        0       -5       0.1        0       10       205       World CoupledTransportation
G4WT0 >     5        0        0      -50       0.1        0       45       250       World biasLimiter
G4WT0 > 1 / 6 DistanceToApply: insideBox -> outsideBox Current Volume: outsideBox Pos: (0,0,-50)
G4WT0 > 1 / 6 GenerateBiasingFinalState: insideBox -> outsideBox Current Volume: parallelWorld Pos: (0,0,-50) -> (0,0,-100)
G4WT0 > 1 from id: 3  to id: 2 outsideBox -> parallelWorld
G4WT0 > 1 Russian roulette: 0.5
G4WT0 >     6        0        0     -100       0.1        0       50       300       World biasLimiter
G4WT0 > 1 / 7 DistanceToApply: outsideBox -> parallelWorld Current Volume: parallelWorld Pos: (0,0,-100)
G4WT0 >     7        0        0     -250       0.1        0      150       450  OutOfWorld CoupledTransportation

To summarize in my view:

At z=50, DistanceToApply is called
Step 3 is calculated, but not limited by biasingLimiter. Instead, boundary of the mass world limits the step.
GenerateBiasingFinalState is called for step 3. PostStepPoint-Position is z=5 → not a biasing boundary, hence no action (“not doing anything”)
Step 4 is calculated, doing physics/transportation in mass world. No biasing world involved here
Step 5 ends at the biasing volume boundary, z=-50
At z=-50, as part of Step 6, DistanceToApply is called again, because last boundary was volume attached to biasing.
Step 6 is performed, moving particle to z=-100 (boundary of outerBox) → Russian roulette

It somehow looks like number 3 in this list prevents from GenerateBiasingFinalState to be called again at the actual boundary for the current biasing volume. Is there a way to make the appropriate manager call GenerateBiasingFinalState again if something else limits the step before the biasing volumes do?

weller · July 3, 2021, 9:08am

hm, this seems to be a complicated issue:frowning:

I have found these slides by @japost (https://indico.cern.ch/event/422168/contributions/1903299/subcontributions/170162/attachments/888120/1249383/EventBiasing.pdf) stating: “Caveats World of importance geometry must ‘overlap’ exactly with mass world”
is this only / also valid for generic biasing? From the context it appears to be “not-generic biasing”.

As an alternative, I tried to use geometric biasing as before, but registering both particles:

  G4String parallelName("myParallelBiasingWorld");
  myParallelWorld* pdet = new myParallelWorld(parallelName);
  detector->RegisterParallelWorld(pdet);
  G4GeometrySampler pgsGamma(pdet->GetWorldVolume(),"gamma");
  pgsGamma.SetParallel(true);

  G4GeometrySampler pgsElectron(pdet->GetWorldVolume(),"e-");
  pgsElectron.SetParallel(true);

  G4VModularPhysicsList* physicsList = new myPhysicsList(0);
  physicsList->RegisterPhysics(new G4ImportanceBiasing(&pgsGamma,parallelName)); // biasing (optional)
  physicsList->RegisterPhysics(new G4ParallelWorldPhysics(parallelName)); // biasing (optional)

  physicsList->RegisterPhysics(new G4ImportanceBiasing(&pgsElectron,parallelName)); // biasing (optional)
  physicsList->RegisterPhysics(new G4ParallelWorldPhysics(parallelName)); // biasing (optional)
  runManager->SetUserInitialization(physicsList);

interestingly, biasing then works perfectly, but only for the particle registered first…
I would like to have both electrons and gamma (and ideally also positrons) to be biased

weller · July 6, 2021, 7:02pm

no ideas, anyone?
Any help is greatly appreciated

weller · July 12, 2021, 11:45am

I am still investigating without success
Does anyone know a solution?

weller · July 13, 2021, 11:28am

maybe a different approach: does anyone know who to ask for help?

I have seen presentation slides from @marcverderi ?

bmorgan · July 13, 2021, 11:40am

Sorry for the lack of response @weller, yes, @marcverderi and also @ahoward would be good people to ask. I know they’ve been very busy with LHCC review documents recently.

ahoward · July 14, 2021, 11:03am

hi @weller, this does look like a bug. I will try to discuss with @marcverderi and get back to you. Out of interest what is the reason that you want to bias more than one particle type at the same time?
Adding multiple particle types in importance biasing is in the work plan. So I could accelerate that if it would solve this issue fastest.

weller · July 14, 2021, 11:50am

hi @ahoward, thanks for your reply.

my intention is to simulate dose distributions in x-ray cameras, to engineer shielding and suitable geometry. radiation comes only indirectly from scintillator to electronics which i want to protect, so I need lots of primary gamma.
Biasing makes my work so much more comfortable, as I can simply multiply secondary/higher order particles by a significant factor only in relevant directions/regions. one alternative would be to simulate the scattering characteristics and use this as a source in a second simulation - geometric biasing allows me to run a all-in-one simulation.

going to much higher intensities and/or high energy x-ray applications (few MeV), e-/e+ seem to start to play a significant role for third degree scattering. biasing gamma by factors of 10-100 per step would skip the electrons, and create events with much higher weight.

I am aware that splitting factors have to be chosen with care, and double check some results with more moderate or no splitting.

weller · July 19, 2021, 8:33am

fyi: I tried running in SerialOnly mode, but same phenomenon.

Is there anything I can do? I could run tests or try things if guided into specific directions… me looking into the code was only marginally useful, the geant4 code is indeed somewhat complex and elaborate

weller · August 3, 2021, 11:44am

hi, just checking back on this… did you, @ahoward, have the change to discuss with @marcverderi?
thanks

weller · September 14, 2021, 6:08am

is there anything new on this issue?

contributing with new content from my side: version 11 has the same issue, but you probably know this anyways

weller · October 13, 2022, 10:54am

hi again
just had a go again with versions up to 11.1-beta1 – same behaviour…

froz1233 · July 26, 2024, 7:47pm

Hello @weller, I wanted to know if you solved your issue.
For information, I have done my own example based on GB03 transforming the parallelworld from slices to cubes (I cut the geometry in the 3 dimensions).
It worked well and I might help you if you still have some issues. I have found a way to bypass the number of particles to bias.

However, when I add a scoring that is mesh-based, the biasing becomes undefined, like if there is a kind of confusion between scoring mesh and biasing mesh. Hence, I checked if the biasing classes have the correct parallel world which is the case. But when I use /tracking/verbose 2, I see that the processName is the name of the scoring mesh and not “biasLimiter” anymore. I am investigating this.