Simulating particle energies beyond 100 TeV

Physics Lists
1

Dear Users,

I am currently trying to simulate particles energies beyond 100 TeV for space based application.
Geant4 version Name: geant4-11-01-beta-01 [MT] (30-June-2022).
I use QGSP_BERT as the physics list.

 G4double minEnergy = 1 * GeV; 
 G4double maxEnergy = 200 * TeV; 
 G4double energy = minEnergy + G4UniformRand() * (maxEnergy - minEnergy);

Unfortunately, when I use particleGun->SetParticleEnergy(energy);

I face the following issues as shown below:

### Run 0 start
--> Event 0 starts.
No model found out of 3
   0.  Elow= 12000, Ehigh= 1e+08   QGSP
   1.  Elow= 3000, Ehigh= 25000   FTFP
   2.  Elow= 0, Ehigh= 6000   BertiniCascade

-------- EEEE ------- G4Exception-START -------- EEEE -------
*** G4Exception : had005
      issued by : G4HadronicProcess::PostStepDoIt
Target element Ca  Z= 20  A= 40
Unrecoverable error in the method ChooseHadronicInteraction of protonInelastic
TrackID= 1  ParentID= 0  proton
Ekin(GeV)= 120615;  direction= (0.794223,0.024219,-0.607144)
Position(mm)= (7.20376e+08,1.41719e+08,-1.49706e+09);  material LunarSurf
PhysicalVolume  <lunar_CrustPhys>
 No HadronicInteraction found out

*** Fatal Exception *** core dump ***
G4Track (0x351376a70) - track ID = 1, parent ID = 0
 Particle type : proton - creator process : not available
 Kinetic energy : 120.615 TeV - Momentum direction : (0.794223,0.024219,-0.607144)
 Step length : 3.39071 cm  - total energy deposit : 17.1894 MeV
 Pre-step point : (7.20376e+08,1.41719e+08,-1.49706e+09) - Physical volume : lunar_CrustPhys (LunarSurf)
 - defined by : hIoni - step status : 4
 Post-step point : (7.20376e+08,1.41719e+08,-1.49706e+09) - Physical volume : lunar_CrustPhys (LunarSurf)
 - defined by : protonInelastic - step status : 4
 *** Note: Step information might not be properly updated.
-------- EEEE -------- G4Exception-END --------- EEEE -------


*** G4Exception: Aborting execution ***

Could you kindly help me out to solve this issue?

Thank you in advance.

1 Like
2

Dear User,
I think that the problem here is that you are setting the energy of the primary particle to a value that lies outside the range covered by the Processes and Models implemented by your Physics List. In fact, if you check the documentation of the QGSP_BERT List, you can see that all the models for the hadronic interactions have an upper limit which, at best, is 100 TeV, while in the event that you are quoting the primary particle features an energy significantly larger (120.615 TeV). This means that Geant4 does not “know” how to determine which processes will occur while tracking the movement of your particle.

The obvious solution is then to correspondingly adjust the upper limits of the Models of all the relevant hadronic and electromagnetic Processes for your physics case. This can be done essentially in two ways:

  1. Harder solution: build your own Physics List, declaring all the relevant Processes and Models and adjusting the upper energy limits with commands such as the following:

EmModel()->SetLowEnergyLimit(minEnergy);
EmModel()->SetHighEnergyLimit(maxEnergy);

  1. Easier solution: using UI commands in your macro, such as (see e.g. here for reference)

/process/eLoss/maxKinEnergy 1 PeV
/process/had/maxEnergy 1 PeV

Honestly I would not recommend so much the first solution, since you could miss a process or two while constructing the List. Using the ReplacePhysics() command could help you not rewriting the entire List, but just the some parts. However, I still think that the second solution is the easier to implement.

One final caveat: pay attention to the fact that merely extending the energy range of your Models may lead to inaccurate results (or, in the worst case scenario, even unphysical results). Check in the literature to know if that is the case.

Best wishes,
Pietro