Proton stuck in loop of cerenkov process with zero step size

Physics Processes, Models and Cross Sections
1

Using 10.7.p1. it was discovered that 6% jobs fail due to a proton getting stuck in a loop of zero step size with cerenkov process involved. This was never observed with the same geometry/physics list etc with previous versions.

Three lines, “PreStep, PostStep, neutronInelastic” belong to one printout from SteppingAction. The numbers are position (x,y,z), time | momentum (x,y,z), Etot | particle name and volume. The third line names the creation process = neutronInelastic, the process of the step, the step length, track length, energy deposit, kinetic energy and track id.

First occurance of track ID: 89651:
PreStep: 1300.6 331.8 238343.4 795.3 | 215.881977 -17.159912 1006.054583 1392.623311 | proton LAV_PbGlBlock
PostStep: 1301.6 331.8 238348.2 795.3 | 214.060463 -22.882968 999.668497 1387.817346 | proton LAV_PbGlBlock
neutronInelastic Cerenkov 4.836132 4.836132 4.805964 449.545333 track ID: 89651

… after a while:

PreStep: 1344.0 310.9 238723.5 797.2 | 70.026040 -5.359012 599.759933 1115.795104 | proton LAV_PbGlBlock
PostStep: 1345.0 310.8 238732.4 797.3 | 42.966892 -7.942397 573.598805 1100.581299 | proton LAV_PbGlBlock
neutronInelastic Cerenkov 8.975197 392.536527 15.213805 162.309286 track ID: 89651
PreStep: 1345.0 310.8 238732.4 797.3 | 42.966892 -7.942397 573.598805 1100.581299 | proton LAV_PbGlBlock
PostStep: 1345.0 310.8 238732.5 797.3 | 45.001214 -7.508665 573.278511 1100.492667 | proton LAV_PbGlBlock
neutronInelastic Cerenkov 0.058150 392.594678 0.088631 162.220654 track ID: 89651
PreStep: 1345.0 310.8 238732.5 797.3 | 45.001214 -7.508665 573.278511 1100.492667 | proton LAV_PbGlBlock
PostStep: 1345.0 310.8 238732.5 797.3 | 44.859313 -7.911088 573.263532 1100.481890 | proton LAV_PbGlBlock
neutronInelastic Cerenkov 0.007718 392.602395 0.010777 162.209877 track ID: 89651
PreStep: 1345.0 310.8 238732.5 797.3 | 44.859313 -7.911088 573.263532 1100.481890 | proton LAV_PbGlBlock
PostStep: 1345.0 310.8 238732.5 797.3 | 45.054419 -7.869496 573.246208 1100.480538 | proton LAV_PbGlBlock
neutronInelastic Cerenkov 0.001587 392.603982 0.001352 162.208525 track ID: 89651
PreStep: 1345.0 310.8 238732.5 797.3 | 45.054419 -7.869496 573.246208 1100.480538 | proton LAV_PbGlBlock
PostStep: 1345.0 310.8 238732.5 797.3 | 44.864572 -7.610301 573.262442 1100.479416 | proton LAV_PbGlBlock
neutronInelastic Cerenkov 0.000817 392.604799 0.001122 162.207403 track ID: 89651
PreStep: 1345.0 310.8 238732.5 797.3 | 44.864572 -7.610301 573.262442 1100.479416 | proton LAV_PbGlBlock
PostStep: 1345.0 310.8 238732.5 797.3 | 44.885146 -7.590507 573.261094 1100.479416 | proton LAV_PbGlBlock
neutronInelastic Cerenkov 0.000179 392.604978 0.000000 162.207403 track ID: 89651
PreStep: 1345.0 310.8 238732.5 797.3 | 44.885146 -7.590507 573.261094 1100.479416 | proton LAV_PbGlBlock
PostStep: 1345.0 310.8 238732.5 797.3 | 44.680343 -7.313612 573.280692 1100.479416 | proton LAV_PbGlBlock
neutronInelastic Cerenkov 0.000179 392.605156 0.000000 162.207403 track ID: 89651
PreStep: 1345.0 310.8 238732.5 797.3 | 44.680343 -7.313612 573.280692 1100.479416 | proton LAV_PbGlBlock
PostStep: 1345.0 310.8 238732.5 797.3 | 44.680338 -7.313611 573.280624 1100.479380 | proton LAV_PbGlBlock
neutronInelastic Cerenkov 0.000179 392.605335 0.000036 162.207367 track ID: 89651
PreStep: 1345.0 310.8 238732.5 797.3 | 44.680338 -7.313611 573.280624 1100.479380 | proton LAV_PbGlBlock
PostStep: 1345.0 310.8 238732.5 797.3 | 44.680304 -7.313606 573.280190 1100.479152 | proton LAV_PbGlBlock
neutronInelastic Cerenkov 0.000158 392.605493 0.000227 162.207139 track ID: 89651
PreStep: 1345.0 310.8 238732.5 797.3 | 44.680304 -7.313606 573.280190 1100.479152 | proton LAV_PbGlBlock
PostStep: 1345.0 310.8 238732.5 797.3 | 44.680299 -7.313605 573.280124 1100.479118 | proton LAV_PbGlBlock
neutronInelastic Cerenkov 0.000029 392.605522 0.000035 162.207105 track ID: 89651
PreStep: 1345.0 310.8 238732.5 797.3 | 44.680299 -7.313605 573.280124 1100.479118 | proton LAV_PbGlBlock
PostStep: 1345.0 310.8 238732.5 797.3 | 44.680299 -7.313605 573.280124 1100.479118 | proton LAV_PbGlBlock
neutronInelastic Cerenkov 0.000009 392.605531 0.000000 162.207105 track ID: 89651
PreStep: 1345.0 310.8 238732.5 797.3 | 44.680299 -7.313605 573.280124 1100.479118 | proton LAV_PbGlBlock
PostStep: 1345.0 310.8 238732.5 797.3 | 44.666283 -7.258610 573.281915 1100.479118 | proton LAV_PbGlBlock
neutronInelastic Cerenkov 0.000009 392.605541 0.000000 162.207105 track ID: 89651
PreStep: 1345.0 310.8 238732.5 797.3 | 44.666283 -7.258610 573.281915 1100.479118 | proton LAV_PbGlBlock
PostStep: 1345.0 310.8 238732.5 797.3 | 44.666283 -7.258610 573.281915 1100.479118 | proton LAV_PbGlBlock
neutronInelastic Cerenkov 0.000009 392.605550 0.000000 162.207105 track ID: 89651
PreStep: 1345.0 310.8 238732.5 797.3 | 44.666283 -7.258610 573.281915 1100.479118 | proton LAV_PbGlBlock
PostStep: 1345.0 310.8 238732.5 797.3 | 44.666283 -7.258610 573.281915 1100.479118 | proton LAV_PbGlBlock
neutronInelastic Cerenkov 0.000009 392.605559 0.000000 162.207105 track ID: 89651
PreStep: 1345.0 310.8 238732.5 797.3 | 44.666283 -7.258610 573.281915 1100.479118 | proton LAV_PbGlBlock
PostStep: 1345.0 310.8 238732.5 797.3 | 44.666283 -7.258610 573.281915 1100.479118 | proton LAV_PbGlBlock
neutronInelastic Cerenkov 0.000009 392.605568 0.000000 162.207105 track ID: 89651
PreStep: 1345.0 310.8 238732.5 797.3 | 44.666283 -7.258610 573.281915 1100.479118 | proton LAV_PbGlBlock
PostStep: 1345.0 310.8 238732.5 797.3 | 44.649307 -7.282236 573.282938 1100.479118 | proton LAV_PbGlBlock
neutronInelastic Cerenkov 0.000009 392.605577 0.000000 162.207105 track ID: 89651
PreStep: 1345.0 310.8 238732.5 797.3 | 44.649307 -7.282236 573.282938 1100.479118 | proton LAV_PbGlBlock
PostStep: 1345.0 310.8 238732.5 797.3 | 44.649307 -7.282236 573.282938 1100.479118 | proton LAV_PbGlBlock
neutronInelastic Cerenkov 0.000009 392.605587 0.000000 162.207105 track ID: 89651
PreStep: 1345.0 310.8 238732.5 797.3 | 44.649307 -7.282236 573.282938 1100.479118 | proton LAV_PbGlBlock
PostStep: 1345.0 310.8 238732.5 797.3 | 44.649307 -7.282236 573.282938 1100.479118 | proton LAV_PbGlBlock
neutronInelastic Cerenkov 0.000009 392.605596 0.000000 162.207105 track ID: 89651
PreStep: 1345.0 310.8 238732.5 797.3 | 44.649307 -7.282236 573.282938 1100.479118 | proton LAV_PbGlBlock
PostStep: 1345.0 310.8 238732.5 797.3 | 44.649306 -7.282235 573.282914 1100.479105 | proton LAV_PbGlBlock
neutronInelastic Cerenkov 0.000009 392.605605 0.000012 162.207092 track ID: 89651
PreStep: 1345.0 310.8 238732.5 797.3 | 44.649306 -7.282235 573.282914 1100.479105 | proton LAV_PbGlBlock
PostStep: 1345.0 310.8 238732.5 797.3 | 44.649305 -7.282235 573.282907 1100.479102 | proton LAV_PbGlBlock
neutronInelastic Cerenkov 0.000002 392.605607 0.000004 162.207089 track ID: 89651
PreStep: 1345.0 310.8 238732.5 797.3 | 44.649305 -7.282235 573.282907 1100.479102 | proton LAV_PbGlBlock
PostStep: 1345.0 310.8 238732.5 797.3 | 44.649305 -7.282235 573.282907 1100.479102 | proton LAV_PbGlBlock
neutronInelastic Cerenkov 0.000000 392.605607 0.000000 162.207089 track ID: 89651
PreStep: 1345.0 310.8 238732.5 797.3 | 44.649305 -7.282235 573.282907 1100.479102 | proton LAV_PbGlBlock
PostStep: 1345.0 310.8 238732.5 797.3 | 44.649305 -7.282235 573.282907 1100.479102 | proton LAV_PbGlBlock
neutronInelastic Cerenkov 0.000000 392.605607 0.000000 162.207089 track ID: 89651
PreStep: 1345.0 310.8 238732.5 797.3 | 44.649305 -7.282235 573.282907 1100.479102 | proton LAV_PbGlBlock
PostStep: 1345.0 310.8 238732.5 797.3 | 44.649305 -7.282235 573.282907 1100.479102 | proton LAV_PbGlBlock
neutronInelastic Cerenkov 0.000000 392.605607 0.000000 162.207089 track ID: 89651
PreStep: 1345.0 310.8 238732.5 797.3 | 44.649305 -7.282235 573.282907 1100.479102 | proton LAV_PbGlBlock
PostStep: 1345.0 310.8 238732.5 797.3 | 44.649305 -7.282235 573.282907 1100.479102 | proton LAV_PbGlBlock
neutronInelastic Cerenkov 0.000000 392.605607 0.000000 162.207089 track ID: 89651

The question is: how to deal with this issue most efficiently. For that it would be great to know what was changed in 10.7.p1 that we now get this issue. Would it be best to take action in the Stepping or StackingAction to kill/prevent such cases? Or already at the user defined physics list with cuts?

I thank you very much in advance for your advice!
Best regards,
Simone

2

Hello,

how for proton tracking neutronsInelastic process appears?
First of all you need to check log output coming at initialisation to be sure that proton and neutron list of processes/models are as desired.

VI

3

Hello,

the physics list was not changed. We just moved to the new G4 version. Hence all processes of interest should be enabled. The initial process is a pion hitting the lead scintillator block, causing a shower. In the relevant steps I ask for the creator process of the proton. It says neutronInelastic, hence I assumed it was created via a neutron interacing inelastically and thus converting to a proton. The process related to the steps is cerenkov.
Simone

4

A few years ago there was this bug:
https://bugzilla-geant4.kek.jp/show_bug.cgi?id=1992
where, when the particle energy was very close to the cutoff energy for Cerenkov radiation, the step length became effectively zero. What’s the index of refraction? 1.9?

Also is the proton energy always the same when the jobs get stuck?

5

The refraction index of the material used is plotted here:
https://refractiveindex.info/?shelf=glass&book=SCHOTT-SF&page=SF57
so, for most wavelengths little below 1.9, i.e. below 1.9 on average, but close to. Yes, the energy is always the same.
Thank you.

6

Could you please share the material composition and density (or Geant4 material name), and refractive index vs. energy that you are using, as ASCII text/code?

Alternatively can you reproduce this with example extended/optical/OpNovice2? So far I haven’t reproduced this.

7

Hello,

I did some investigations. But I am not able to reproduce this issue with a simpyfied version of the detector: I used the same geometry of the detector part concerned, the same material and the same initial proton (BWT: the proton is produced by a pi- hitting another detector farer away, there, producing a neutron that finally creates the proton via neutroninelastic in the lead glass. However, the parent ID of the proton is that of the pion, which I find misleading since it was directly produced by the neutron converting to a proton.)

The point is: we never saw this issue before moving from 10.6.p2 to 10.7.p1. We changed neither the material nor the detector geometry of the detector concerned. Hence the issue must be related to a change in the 10.7.p1. Either it was not spotted on our side with the older version since it was “hidden” and is now revealed by a change in g4 or it is a real problem in the new version. In any case we have to understand what is going on, otherwise we cannot use the latest version.

Do you have any idea what could be checked further? I.e. which printouts?

One thing I observed in the real detector is in SteppingAction.cc: aStep->GetNumberOfSecondariesInCurrentStep() is always 0, while asking G4Cerenkov class directly before finalising the process the step indeed has secondary photons produced (from poissonian distribution). This is not the case for the toy example.
Furthermore: after step limit of 0.00000000000001 is reached, the proton is ‘stuck’, the min, ekin = 162.207 and beta = 0.522559. In the toy example, however, the Cerenkov process is stopped at smaller step limits, with finally 0 photons and same ekin,beta but ionisation or protoninelastic takes place right afterwards. This does not happen in the full detector simulation.

All details of the lead glass can be found in the above link but I can also send you the material as it is built in G4 code convention if you prefer.

Thanks a lot for any help regarding this issue in advance.

8

Are you set up to re-compile the Geant4 source code? If so could you please try changing this line:

https://geant4.kek.jp/lxr/source/processes/electromagnetic/xrays/src/G4Cerenkov.cc#L453

if(Step < 1.e-15 * mm)
  return StepLimit;

so that the value in the condition is greater than the step size you are seeing when the particle is stuck.

9

Thank you for the hint. If I change the condition to <=1e-14mm as in the case I observe above, then indeed everything works well, i.e. the stepping is progressing and terminating properly.
I wonder why this happens, since I would guess there are other mechanisms to kill the process/particle if not advancing, right?
The other question from my previous reply is: why do I see 0 secondaries when asking for this info in the SteppingAction even for larger step lengths while from G4Cerenkov the number of photons is indeed larger 0 at the same time. This happens only for the real detector design. It is the case for all Cerenkov steps of the proton, so also for the very first (large step lengths) steps.

10

Concerning the last point: Do you have any idea what else I could check? To which class the is the number of photons from G4Cerenkov handed over before it is accessible via aStep->GetNumberOfSecondariesInCurrentStep() in the SteppingAction?

11

In the OpNovice2 example there is

const std::vector<const G4Track*>* secondaries =
      step->GetSecondaryInCurrentStep();

which seems to work. Maybe add this to your code, get the length of the vector, and compare to

aStep->GetNumberOfSecondariesInCurrentStep()

They should be the same.

12

Thank you for the hint.
Yes, they are the same, both yielding 0 although G4Cerenkov has photons registered. So we have three issues here:

  1. the stuck proton: changing the steplimit in G4Cerenkov is not an option I think, since we would like to use standard G4 version, hence the issue would need to get understood better: how, in which class, which parameters?
  2. the parentID of the proton being that of the first parent, the pion, although it should be that of direct mother, which is the neutron that is created by the pion
  3. the number of secondaries being 0 although there some photons being produced.

Concerning 3): wherelse could I place a printout to see where the info of photons gets lost? In the process manager?

Concerning 1): In the meantime I also run with 10.7.p2, the same, 6% of the jobs are killed due to wall time exceeded since the proton gets stuck.

Thank you very much.

13

I don’t know what is causing this. I suspect something numerical (i.e. a very small number ‘equals’ zero). I could adjust the parameter (again) but that’s not very satisfying. Or expose the parameter to the user. Ideally I’d understand the problem better. This is clearly a bug–if you like, please open a bug report.

This is a hadronics question, and I don’t know the answer. I suggest starting a new discussion specific to this.

I don’t know. Do you call that method in your UserSteppingAction? Are you sure you’re checking the correct particle? Maybe try
/tracking/verbose 2 (or 3?)
to list the secondaries for each step.

14

at 3) from above is understood: a complex chain of selections/cuts finally leads to killing optical photons to advance further in the concerned geometry from our side. But this does not seem to have an effect on 1), since this part in our code was not touched and I tested it by switching this on/off.

at 1) The question is still, what could have changed in 10.7. such hat the standard killing mechanisms in G4 in case the mother proton is not advancing and not producing photons any longer (checked with G4Cerenkov printout) is not always effective. Ok, I will finally open a bug report and I will also mention 2) there. We had some delay due other important issues to be solved and vacation period, hence the late reply.
Thank you!

15

I think the answer to 1) is to set the value here:
https://geant4.kek.jp/lxr/source/processes/electromagnetic/xrays/src/G4Cerenkov.cc#L453

to

  if(Step < G4ThreeVector::getTolerance())
    return StepLimit;

(which is about 2.2e-14 mm).

Unless there is feedback that this doesn’t work, I’ll change it for version 11 and 10.7.patch3 (if there will be one).

16

Yes, this would help. Great, thank you!