HepMC3 reader with Pythia8 Argantyr generator: status checker not implemented

[espadaro1]

We are trying to use an input file in the HepMC3 format. This file is generated with Pythia8 Argantyr. When we run it with ddsim, we got the following error:

G4DecayTable::SelectADecayChannel :: no possible DecayChannel      neutron

-------- EEEE ------- G4Exception-START -------- EEEE -------
*** G4Exception : DECAY003
     issued by : G4Decay::DoIt
Can not determine decay channel for neutron
 mass of dynamic particle: 0.93827 (GEV)
 dacay table has 1 entries
0: BR 1, IsOK? 0, --> e- + anti_nu_e + proton

*** Fatal Exception *** core dump ***
G4Track (0x126bfb28) - track ID = 38, parent ID = 0
Particle type : neutron - creator process : not available
Kinetic energy : 0 eV  - Momentum direction : (0,0,-1)
Step length : 0 fm  - total energy deposit : 0 eV
Pre-step point : (0,0,0) - Physical volume : BeamPipeVacuum_62 (Vacuum)
 - defined by : not available
 Post-step point : (0,0,0) - Physical volume : BeamPipeVacuum_62 (Vacuum)
- defined by : not available
*** Note: Step information might not be properly updated.
-------- EEEE -------- G4Exception-END --------- EEEE -------

*** G4Exception: Aborting execution ***
 *** Break *** abort

It seems that the hepmc3 reader does not check for the particles status codes and passes all particles to Geant4, also the ones with status code 3 (documentation) and 4 (beam particles). Indeed, the neutron does not have the correct mass because it is a virtual particle in the Pythia8 event record. We have seen that in the HepMC2 reader a status code checker was implemented.

Here, I give the details to reproduce the error:

• lxplus centos7, GCC 11
• DD4hep version: v1.20
• Run ddsim with SiD.xml compact file, as an example:

  ddsim --dumpSteeringFile > steering.py
  ddsim --steeringFile steering.py --compactFile DDDetectors/compact/SiD.xml --inputFiles Pythia_output.hepmc --outputFile MC_Py_10.root --hepmc3.useHepMC3 True -N 10 --part.userParticleHandler='' -vDEBUG > output.log

I attach the input file, where we removed particles with status >10 and the output log. output.log Pythia_output.hepmc.txt --> Rename to Pythia_output.hepmc

Thank you for your support

[andresailer]

Please have a look at #920

[andresailer]

I think we have a problem with the nuclei (1000741840 and 1000711809 in this file), if I change them to protons the crash disappears.

EDIT: will investigate to enable the necessary Physics parts

[andresailer]

Hi @espadaro1 I now have a working simulation running in #920 One issue that showed up was the excited state of the Lutetium Ion, I am still trying to figure out if we can give that to geant4, or have to pretend it isn't excited as I am doing at the moment.

https://github.com/AIDASoft/DD4hep/pull/920/files#diff-2dcb1a8b3bf461264048963c9ff6e142d8334e9389e9567e9062717e4b581fa8R123-R124

What do you think?

[andresailer]

One important commit I added only 2 hours ago: https://github.com/AIDASoft/DD4hep/pull/920/commits/26683ca4120c5311597b620daf2c7a34666e1052

[andresailer]

I tried a different G4Ion constructor

    int id = particle->pdgID;  // is 10 L ZZZ AAA I
    int L, A, Z, I;
    id -= 1000000000;  // leading 10 is just identifier
    L = id / 10000000; id %= 10000000;
    Z = id / 10000;    id %= 10000;
    A = id / 10;       id %= 10;
    I = id;
    return G4IonTable::GetIonTable()->GetIon(Z, A, L, I);

But that also gives

-------- WWWW ------- G4Exception-START -------- WWWW -------
*** G4Exception : PART105
      issued by : G4IonTable::GetIon()
Ion cannot be created by an isomer level. Use excitation energy.
*** This is just a warning message. ***
-------- WWWW -------- G4Exception-END --------- WWWW -------

So I don't think there is sufficient information in the MC record to use the excited ion in the simulation?

[espadaro1]

Thanks Andre for checking this. Yes, probably there is some missing information in the MC record for excited ions. Do you know which information is needed to create ions in Geant4? We can find out if this information is available in Pythia.

Apart from the excited ions, it seems that the output is correct. We will also try to process more events to see if there are no other errors.

[andresailer]

There are some "GetIon" functions taking the G4Double E argument https://geant4.kek.jp/Reference/v11.0.2/classG4IonTable.html

[MarkusFrankATcernch]

Is there any further update needed for this issue? If not, please close the issue.

[espadaro1]

Hi Markus,

sorry, this issue is open because we are still trying to understand how to treat excited ions and the solution proposed is only partial and has not been merged into the main project.

About the issue, it seems strange that there is no possibility to compute the excited energy from the excitation level in Geant. No particular information should be passed from the generator file.

These ions have an excited level which is set as ambiguous by Geant4. Maybe this is why it cannot recover the energy. Do you have suggestions about who Geant experts I could contact for this issue?

Thank you.

Best regards, Elisabetta Spadaro

Il giorno 27 ott 2022, alle ore 15:37, MarkusFrankATcernch @.**@.>> ha scritto:

Is there any further update needed for this issue? If not, please close the issue.

— Reply to this email directly, view it on GitHubhttps://github.com/AIDASoft/DD4hep/issues/918#issuecomment-1293539228, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AZVAKVVWBPRPUXQE2HZRMMDWFKAQZANCNFSM5Y6APMLQ. You are receiving this because you were mentioned.Message ID: @.***>

[andresailer]

Probably then a question for https://geant4-forum.web.cern.ch/

[espadaro1]

I think I have solved the issue, before asking to experts. Sorry for taking so long.

This is the proposed modification inside the method Geant4ParticleHandle::definition() of Geant4Particle.cpp:

const G4ParticleDefinition* Geant4ParticleHandle::definition() const   {
  G4ParticleTable*      tab = G4ParticleTable::GetParticleTable();
  G4ParticleDefinition* def = tab->FindParticle(particle->pdgID);
  if( 1000000000 < particle->pdgID) {
    int id = particle->pdgID;  // is 10 L ZZZ AAA I
    int L, A, Z, lvl;
    id -= 1000000000;  // leading 10 is just identifier
    L = id / 10000000; id %= 10000000;
    Z = id / 10000;    id %= 10000;
    A = id / 10;       id %= 10;
    lvl = id;

    G4double E = 0.0;
    G4IonTable*      tab_ion = G4IonTable::GetIonTable();
    G4ParticleDefinition* def_ion = tab_ion->FindIon(Z, A, lvl);

    if(lvl>0){
      if(def_ion != NULL ) {
        E = ((const G4Ions*)(def_ion))->GetExcitationEnergy();
        //std::cout<<"IsomerLevel, "<<lvl<<", excitation energy E="<<E<<std::endl;
      }
    }
    return G4IonTable::GetIonTable()->GetIon(Z,A,L,E,0);
  }
  if ( 0 == def && 0 == particle->pdgID )   {
    if ( 0 == particle->charge )
      return G4Geantino::Definition();
    return G4ChargedGeantino::Definition();
  }
  return def;
}

Let me know if it is works for you as well and if you need anything else from me.

[andresailer]

Hi @espadaro1 Could you open a PR with all the changes to make it work for you? I am not sure on top of which developments to apply your proposed changes.

Thanks, Andre

[elspadar]

Hi @andresailer, the changes can be applied on top of the pull request #920. I was working on your rejectGen branch. There, you modified Geant4Particle.cpp, and what I proposed can go there.

In addition, I have just noticed that in the previous PR you also modified the class Geant4ParticleHandler.cpp, setting the last digit of the pdgID of ions to 0, in line 347. I am testing if it can work without that modification now. Therefore, I am changing line 347 to:const int pdgID = p->pdgID;. However, it is taking too long time to process the events, so I wonder why you included that change in the previous version and what happens if we keep the previous implementation. We are not associated the correct G4PrimaryParticle for the excited ions, right? If you can help on this final issue, I would appreciate. Thank you.

So I think that, if you can, is better if you apply these changes directly in your rejectGen branch.

[andresailer]

Hi @elspadar , The change for the pdgID was needed to align the particle information when ignoring the excitation level, otherwise the program crashed in the end (I think, it was too long ago). As we can keep the excitation level this change has to go out again.

But I also run into problems when simulating with the changes shown on top, so I wanted to be sure from what starting point to go.

[andresailer]

Hi @elspadar @espadaro1 Does this work for you?

    G4ParticleDefinition* def_ion = tab_ion->FindIon(Z, A, lvl);

I don't see this ever returning an ion when lvl > 0. I only see nullptrs

[elspadar]

Hi @andresailer,

yes, it works for me. I can get ions with lvl>0 and with associated excitation energy. The class G4IonTable has been updated in latest versions of Geant4 (I think from 10.3.2 but I cannot find it right now), with new implementations of FindIon and GetIon methods. Could it be a reason of nullptrs? I am using Geant4-10.7.2 right now.

[andresailer]

I am also using Geant4 version Name: geant4-10-07-patch-02 [MT] (11-June-2021)

Which Physics List are you using? I am running with FTFP_BERT

[elspadar]

Ok. I am also using FTFP_BERT. I have also tested different physics lists and I can run with everyone.

I remembered I had similar problems in the past with null pointers but with the code I provided above it worked and I do not think I have changed anything else in dd4hep. But I will try think about it again.

[andresailer]

If you could push dd4hep as you have it, we can directly compare the code.

[andresailer]

So the excitation level works in some cases (I added some printouts to get pdgID and pointer(or nullptr)

Creating ION 1000731763
Found Ion ? 0x15a42740
IsomerLevel, 3, excitation energy E=2.771

But not for the ion from the Monte Carlo record

Creating ION 1000711809
Found Ion ? 0
Returning ion with PDG 1000711800

And then the excitation level is zeroed similar to what I was doing before. Is it working for your for this ion 1000711809 ?

[espadaro1]

Ok. For this ion, I also find a null pointer. So I think we get the same results. It seems there is no excitation energy associated to some ions. Maybe this is something to ask on the geant4 forum?

[andresailer]

Yes, I think we need their help in this case.

[espadaro1]

Ok, I have submitted a new request on the geant4 forum: https://geant4-forum.web.cern.ch/t/g4iontable-findion-with-excitation-level-nullptrs/9224. Let's see what they reply. Thanks.