Geant4InputHandling: Treat Ions with Excitation Energy?

DDG4/python/DDSim/Helper/Physics.py

@@ -5,8 +5,11 @@
 
 from DDSim.Helper.ConfigHelper import ConfigHelper
 from g4units import mm
+import logging
 import ddsix as six
 
+logger = logging.getLogger(__name__)
+
 
 class Physics(ConfigHelper):
   """Configuration for the PhysicsList"""
@@ -139,7 +142,11 @@ def setupPhysics(self, kernel, name=None):
       rg.RangeCut = self.rangecut
 
     for func in self._userFunctions:
-      func(kernel)
+      try:
+        func(kernel)
+      except Exception as e:
+        logger.error("Exception in UserFunction: %r", e)
+        raise RuntimeError("Exception in UserFunction: %r" % e)
 
     return seq
 

DDG4/src/Geant4InputHandling.cpp

@@ -21,10 +21,10 @@
 #include <CLHEP/Units/PhysicalConstants.h>
 
 // Geant4 include files
-#include <G4ParticleDefinition.hh>
 #include <G4Event.hh>
-#include <G4PrimaryVertex.hh>
+#include <G4ParticleDefinition.hh>
 #include <G4PrimaryParticle.hh>
+#include <G4PrimaryVertex.hh>
 
 // C/C++ include files
 #include <stdexcept>
@@ -344,7 +344,9 @@ int dd4hep::sim::smearInteraction(const Geant4Action* caller,
 static G4PrimaryParticle* createG4Primary(const Geant4ParticleHandle p)  {
   G4PrimaryParticle* g4 = 0;
   if ( 0 != p->pdgID )   {
-    g4 = new G4PrimaryParticle(p->pdgID, p->psx, p->psy, p->psz, p.energy());
+    // For ions we use the pdgID of the definition, in case we had to zero the excitation level, see Geant4Particle.cpp
+    const int pdgID =  p->pdgID < 1000000000 ? p->pdgID : p.definition()->GetPDGEncoding();
+    g4 = new G4PrimaryParticle(pdgID, p->psx, p->psy, p->psz, p.energy());
   }
   else   {
     const G4ParticleDefinition* def = p.definition();

DDG4/src/Geant4Particle.cpp

@@ -12,17 +12,20 @@
 //==========================================================================
 
 // Framework include files
-#include <DD4hep/Printout.h>
-#include <DD4hep/Primitives.h>
 #include <DD4hep/InstanceCount.h>
+#include <DD4hep/Primitives.h>
+#include <DD4hep/Printout.h>
 #include <DDG4/Geant4Particle.h>
-#include <TDatabasePDG.h>
-#include <TParticlePDG.h>
-#include <G4ParticleTable.hh>
-#include <G4ParticleDefinition.hh>
-#include <G4VProcess.hh>
+
 #include <G4ChargedGeantino.hh>
 #include <G4Geantino.hh>
+#include <G4IonTable.hh>
+#include <G4ParticleDefinition.hh>
+#include <G4ParticleTable.hh>
+#include <G4VProcess.hh>
+
+#include <TDatabasePDG.h>
+#include <TParticlePDG.h>
 
 #include <sstream>
 #include <iostream>
@@ -114,6 +117,32 @@ void Geant4Particle::removeDaughter(int id_daughter)  {
 const G4ParticleDefinition* Geant4ParticleHandle::definition() const   {
   G4ParticleTable*      tab = G4ParticleTable::GetParticleTable();
   G4ParticleDefinition* def = tab->FindParticle(particle->pdgID);
+  if( 1000000000 < particle->pdgID) {
+    // creating ions here
+    // ion encoding is 10 L ZZZ AAA I
+    int id = particle->pdgID - 1000000000;  // leading 10 is just identifier
+    const int L = id / 10000000; id %= 10000000;  // strange-quark content
+    const int Z = id / 10000;    id %= 10000;
+    const int A = id / 10;       id %= 10;
+    const int lvl = id;
+    G4IonTable* tab_ion = G4IonTable::GetIonTable();
+    // Have to look for existing Ions, excited or not
+    G4ParticleDefinition* def_ion = tab_ion->FindIon(Z, A, L, lvl);
+    if(def_ion) {
+      //We found an existing Ion, we are good to go
+      printout(VERBOSE,"Geant4Particle","+++ Returning ion with PDG %10d", def_ion->GetPDGEncoding());
+      return def_ion;
+    } else if(lvl == 0) {
+      // GetIon creates the Ion if it does not exist, if this does not work something is seriously wrong
+      printout(VERBOSE,"Geant4Particle","+++ Creating ion with PDG %10d", particle->pdgID);
+      return tab_ion->GetIon(Z, A, L, 0.0);
+    }
+    //Cannot use GetIon with lvl > 0, must give energy, but we don't know where to get energy from
+    printout(WARNING,"Geant4Particle","+++ Cannot find excited ion with PDG %10d, setting excitation level to zero",
+             particle->pdgID);
+    return tab_ion->GetIon(Z, A, L, /* E= */ 0.0);
+  } // finished with ions
+
   if ( 0 == def && 0 == particle->pdgID )   {
     if ( 0 == particle->charge )
       return G4Geantino::Definition();

DDTest/CMakeLists.txt

@@ -57,6 +57,16 @@ if (DD4HEP_USE_GEANT4)
     pytest ${PROJECT_SOURCE_DIR}/DDTest/python/test_import_ddg4.py)
   SET_TESTS_PROPERTIES( t_test_python_import_ddg4 PROPERTIES FAIL_REGULAR_EXPRESSION  "Exception;EXCEPTION;ERROR;Error" )
 
+  if(DD4HEP_USE_HEPMC3)
+    ADD_TEST( t_test_ddsim_ion_hepmc3 "${CMAKE_INSTALL_PREFIX}/bin/run_test.sh"
+      ddsim --compactFile=${CMAKE_INSTALL_PREFIX}/DDDetectors/compact/SiD.xml --runType=batch -N=1
+      --outputFile=testSidIons.root --inputFiles ${CMAKE_CURRENT_SOURCE_DIR}/inputFiles/Pythia_output.hepmc
+      --part.userParticleHandler=)
+    # Tests causes G4Exception about: Isomer level 9 may be ambiguous, just a warning
+    # space before Exception intentional!!!
+    SET_TESTS_PROPERTIES( t_test_ddsim_ion_hepmc3 PROPERTIES FAIL_REGULAR_EXPRESSION  "ERROR;Error; Exception" )
+  endif()
+
   foreach(TEST_NAME
       test_EventReaders
       )