Physics list for decay chain and gamma cascades

Getting Started
8

I think I may have found the problem now and it seems to be a bug in Geant4’s code. I will do my best to summarize from the beginning what I think happens.

  1. I am simulating nuclear decays where gammas are emitted. It is important for me that the gammas are correlated and that the angular distributions of gamma-ray cascades are accurate. Therefore, in my PhysicsList.cc, I am using:
    G4DeexPrecoParameters* deex = G4NuclearLevelData::GetInstance()->GetParameters();
    deex->SetCorrelatedGamma(true);
  1. Geant4 simulates a transition from an excited nuclear state to a lower one or a ground state. The states each have an angular momentum (J) and parity (P). The data for this is stored in YourInstallFolder\share\Geant4\data\PhotonEvaporation5.7. Here J and P information is stored as J*P. For example, J=2, P=+ is stored as 2 and J=3/2, P=- is stored as -1.5.
  2. The problem arises when the second nuclear state has negative parity. So P=- and the value for J*P is <0. Let us call this JP2 < 0.
  3. Geant4 uses Geant4SourceFolder\source\processes\hadronic\models\de_excitation\photon_evaporation\src\ to handle correlated gammas. Somewhere here JP2 is multiplied by two.
  4. In G4PolarizationTransition.cc, our 2*JP2 is stored in variable fTwoJ2. In the case above, it is negative. With verbosity, it is printed and valued as expected.
  5. This happens on line 279 in G4PolarizationTransition.cc:
  std::size_t newlength = fTwoJ2+1;
  1. But std::size_t is unsigned and fTwoJ2+1 can be negative. In the problematic case newlength caps out at ~2^64. This is seen in frame #0 in my last post where newlength=18446744073709551611.
  2. On line 287 in G4PolarizationTransition.cc there is a for loop over newlength, which now is very big:
for(G4int k2=0; k2<(G4int)newlength; ++k2) {
  1. From here on, bad stuff happens which I think causes the segmentation fault.

I have tested several examples and this segfault happens every time fTwoJ2+1<0. Also fTwoJ2+1=0 causes a segfault. So this segfault always happens when the parity of the second state is negative. It never happens when the parity of the second state is positive or J=0.

There are many nuclear decays where this happens.

It feels strange that this has gone undetected and I am not excluding the possibility that I am wrong. I am trying to learn both Geant4 and C++ here, so please inform me of any mistakes.