Replace part of reference physics list

Physics Lists
1

First apologies if this is answered somewhere else but I have spent some time going through the manuals, examples and forum posts and can not work out the answer to this. I have a model which applies just to inelastic collisions of pions in a limited energy range and is available for certain elements. I want to have a default physics list such as QGSP_BERT for everything else. The model was written by someone else but so far I can see it works but I can’t work out how to overlay it onto a reference physics list. I tried this:
G4PhysListFactory* factory=new G4PhysListFactory();
G4VModularPhysicsList* PhysicsList=factory->GetReferencePhysList(“QGSP_BERT”);
PhysicsList->RegisterPhysics(new MyPionPhysicsModel());
runManager->SetUserInitialization(PhysicsList);

but as soon as the generator produces a particle outside the energy range of the model it crashes saying

G4EnergyRangeManager:GetHadronicInteraction: counter=1, Ek=5861.95, Material = G4_C, Element = C
0 low=0, high=5995
In /t2k/Software/13.2/Geant4_10.1.03.01/geant4_source/source/processes/hadronic/management/src/G4EnergyRangeManager.cc, line 136:
===> GetHadronicInteraction: No Model found

-------- EEEE ------- G4Exception-START -------- EEEE -------
*** G4Exception : had005
issued by : G4HadronicProcess::PostStepDoIt
In /t2k/Software/13.2/Geant4_10.1.03.01/geant4_source/source/processes/hadronic/management/src/G4EnergyRangeManager.cc, line 136:
===> GetHadronicInteraction: No Model found
Target element C Z= 6 A= 12
Unrecoverable error in the method ChooseHadronicInteraction of MyPionPhysicsModel
TrackID= 1 ParentID= 0 pi+
Ekin(GeV)= 5.86195; direction= (2.46853e-05,-1.52839e-05,1)
Position(mm)= (1.18246e-05,-3.72635e-06,-0.582798); material G4_C
PhysicalVolume
No HadronicInteraction found out

*** Fatal Exception *** core dump ***
-------- EEEE -------- G4Exception-END --------- EEEE -------

*** G4Exception: Aborting execution ***

What is the correct way to do this?

2

Hello,

you cannot add a single model via this interface. The exception is correct.

For today, the only possibility is to create your own physics constructor and replace G4HadronPhysicsQGSP_BERT. Such interface to replace a component of a physics list works fine. You may take G4HadronPhysicsQGSP_BERT as a starting point, rename it locally, and inside add your model to the pion process.

VI

3

ok thanks I’ll give that a try.

4

I’m getting some confusing results. As a test I’ve tried creating a physics list using
GetReferencePhysList(“QGSP_BERT”);
and then replacing physics with
PhysicsList->ReplacePhysics(new G4HadronPhysicsQGSP_BERT())
which I exected to be identical, however I get different results when I analyse what is produced.
Is this expected?
In the print out for hadronic proceses I get duplicated sections like this:
Process: pi+Inelastic
Model: QGSP: 12 GeV —> 100 TeV
Model: FTFP: 9.5 GeV —> 25 GeV
Model: BertiniCascade: 0 eV —> 9.9 GeV
Cr_sctns: G4CrossSectionPairGG: 0 eV —> 100 TeV
G4CrossSectionPairGG: G4PiNuclearCrossSection cross sections
below 91 GeV, Glauber-Gribov above
Cr_sctns: G4CrossSectionPairGG: 0 eV —> 100 TeV
G4CrossSectionPairGG: G4PiNuclearCrossSection cross sections
below 91 GeV, Glauber-Gribov above
Cr_sctns: GheishaInelastic: 0 eV —> 100 TeV

Process: pi+Inelastic
Model: QGSP: 12 GeV —> 100 TeV
Model: FTFP: 9.5 GeV —> 25 GeV
Model: BertiniCascade: 0 eV —> 9.9 GeV
Cr_sctns: G4CrossSectionPairGG: 0 eV —> 100 TeV
G4CrossSectionPairGG: G4PiNuclearCrossSection cross sections
below 91 GeV, Glauber-Gribov above
Cr_sctns: G4CrossSectionPairGG: 0 eV —> 100 TeV
G4CrossSectionPairGG: G4PiNuclearCrossSection cross sections
below 91 GeV, Glauber-Gribov above
Cr_sctns: GheishaInelastic: 0 eV —> 100 TeV