added lightIonQMD as a special physics constructor

[GitFuchs]

This is based on a discussion with Nils Krah. The full email discussion is added below.

Summary: In Geant4 11.2.x a new physics option G4LightIonQMDPhysics became available (see Sato et al., PMB 2022(67) DOI:10.1088/1361-6560/ac9a9a). It is supposed to yield better results with fragmentation in the clinical energy range.

It is already included into the physics list Shielding, but due to an incomplete implementation can not yet be used with the default physics factory from Geant4.
A report to the Geant4 developers is already pending https://bugzilla-geant4.kek.jp/show_bug.cgi?id=2615

As a workaround a special physics constructor G4LightIonQMDPhysics was made available.
Invoking it, will unload already registered physics for these processes (e.g. standard QMD) and replace it with G4LightIonQMDPhysics.

A small description was added to the documentation as well.

############### Email exchange below
31th of May
Dear all,
I found out that the new Geant4.11.2.x has a new physics option, G4LightIonQMDPhysics. It is supposedly better when simulating fractionation tails for carbon beam therapy...
It can be used with the physics lists Shielding if an additional option in the physcis list is set.

Shielding::Shielding(G4int verbose,const G4String& n_model,const G4String& HadrPhysVariant, G4bool useLightIonQMD)

Unfortunately, the physics list factory of Geant4 does not integrate it, yet.
Ideally in G4PhysListFactory.cc an additional line should be input such as
else if(had_name == "Shielding_LQMD") {p = new Shielding(verbose,"","",True);}
Such an option may come, but I would like to get it into OpenGATE right away ;-)
What I would like to do is to create a method to call the physics list with options from OpenGATE
In c++ the following lines would be required:
G4VModularPhysicsList* p = nullptr;
p = new Shielding(verbose,"LEND","",True); //the hadronic part
p->ReplacePhysics(new G4EmStandardPhysics_option4(verbose)); //the EM part
The ReplacePhysics part is already in OpenGATE (see pyG4ModularPhysicsList.cpp)
How could I approach the calling of the class (including options)?
In managers.py I found the function

create_modular_physics_list_class
(
g4_physics_constructor_class_name
)
But it does not seem to be doing what I want to achieve.
Could anyone help me out here?
Could anyone help me out here?
All the best,
Hermann

4th of June
Hey Hermann,
​thanks for lookinig into this. I remember you already mentioned this issue during the Gate meeting 2 weeks ago.

​The way to activate this new option seems poorly/incompletely implemented on the Geant4 side because it is not compatible with the PhysicsListFactory. Ideally, G4 should fix this, but for the time being, I guess we (GATE) need to find a reasonable workaround.
​
​Quick explanation to start with:
GATE creates the physics list in PhysicsListManager.get_physics_list() given the name provided by the user. ​
​
​The function create_modular_physics_list_class ​is actually made for PhysicsConstructor classes.
​Reason:
​We want Gate users to be able to use, e.g., G4EmStandardPhysics as a physics list, but technically this is not G4ModularPhysicsList, but a G4PhysicsConstructor. The create_modular_physics_list_class creates a physics list of the same name and adds the physics constructor to it. This homogenizes the code.
​
​The Shielding list is really a physics list and GATE creates it via the G4PhysicsListFactory and there we run into the problem you mention.

​
​Regarding the issue, ​I see multiple alternative workarounds:
​
​1) In PhysicsListManager.get_physics_list(), capture the name "ShieldingLQMD" and call the constructor of the G4 Shielding class directly from python (needs to be wrapped) with the boolean variable set to true.
2) ​Similar to 1, but implement a simple C++ function which creates a pointer to the physics list, calls "p = new Shielding(..., true);" and returns the pointer. In PhysicsListManager.get_physics_list(), call this function.
​3) Implement a GatePhysicsListFactory class which inherits from G4PhysicsListFactory and simply intercepts the input variable. If the variable is "ShieldingLQMD", create the physics list manually with "p = new Shielding(..., true)", if not, call the method from the base class.
​
​Option 3 is similar to option 2, but handles the logic at the C++ level. The python code remains untouched.
​
​My personal suggestion:
​Let's reach out to the G4 devs and ask if they plan to "fix" this issue and when. Then we decide which workaround to add.
​
​If you want to work on this, feel free to drop a line and we can look at the details together.
​
​Best,
​Nils

5th of June
Dear Nils,
I agree with you, this is an implementation issue of Geant4. I have opened a ticket, let's see what the answer is.
Thank you very much for providing your thoughts! Personally, I do like to include something in the official GATE code base to fix an inconsistent programming on the Geant4 side, so I will keep this in a personal dev branch.

I might have found an additional option, which might also be interesting, especially if one wants to tune physics lists. However, I am not sure if it really does what I think it does.
It seems that the special_physics_constructors also provide an option to replace modular physics lists (see managers.py line 416)
So in the class PhysicsListManager
line 332 ff
class PhysicsListManager(GateObject):
available_g4_physics_constructors = [...
I added G4LightIonQMDPhysics
This way, after selecting Shielding as physics list, using
sim.physics_manager.special_physics_constructors.G4LightIonQMDPhysics = True
should replace the default G4IonQMDPhysics of Shielding with G4LightIonQMDPhysics.

Looking at the output, the physics does change, but I am not sure whether this hack is actually doing what it should.

What do you think?
All the best,
Hermann

5th of June
Dear Hermann,
​your option is really good and much more efficient than the ones I proposed.
​
​I had not looked at the (new) Shielding class in detail before, but did so now.
​In Shielding.cc in Geant4 11.2.1, I find:

​​ if (useLightIonQMD){
RegisterPhysics( new G4LightIonQMDPhysics(verbose) );
} else {
RegisterPhysics( new G4IonQMDPhysics(verbose) );
}

​useLightIonQMD is false by default. So either of the two PhysicsConstructors is registered.
​
​The user input parameter special_physics_constructors (a dictionary/Box actually) is, as you correctly identified, an interface for the user to activate additional PhysicsConstructors (including optical, chemistry, etc). Internally, GATE calls the wrapped G4 function ReplacePhysics() which in turn checks the type of PhysicsConstructor. If such a type of physics is already present in the physics list, it is replaced. If not, the physics constructor is added.
​
​In the specific case, G4LightIonQMDPhysics and G4IonQMDPhysics, both have physics type "bIons" (value 8), so G4IonQMDPhysics is really replaced if G4LightIonQMDPhysics is activated as additional constructor via the user parameter described above.
​Therefore, the result of the mechanism is identical to the logic of the lines in Shielding.cc (see above).
​
​​To implement the additional G4LightIonQMDPhysics physics constructor, you need to do 2 things:
​(probably what you already did)
​

1. add a line
   ​ADD_PHYSICS_CONSTRUCTOR(G4LightIonQMDPhysics)
   ​in the wrapper code in pyPhysicsLists.cpp
   ​
2. add an entry to the class-level dictionary in PhysicsListManager (python)
   ​ special_physics_constructor_classes["G4LightIonQMDPhysics"] = g4.G4LightIonQMDPhysics
   ​​The PhysicsManager then automatically picks up this constructor as an option and makes it available to the user via
   ​sim.physics_manager.special_physics_constructors.G4LightIonQMDPhysics = True
   ​as you said.
   ​
   ​Can you create a PR so we can merge the change into master?
   ​
   ​Best,
   Nils​