AIDASoft / DD4hep

Detector Description Toolkit for High Energy Physics
http://dd4hep.cern.ch
GNU Lesser General Public License v3.0
49 stars 95 forks source link

Protect against wrong primary particle mass #1233

Closed atolosadelgado closed 5 months ago

atolosadelgado commented 7 months ago

Hi,

ddsim allows to load primary particles from an external file, but it seems there is no sanity check for the mass of the particles.

Question: would be worth to check the mass of the primary particles before starting the simulation? This issue would affect the builtin particle gun, which allows to define a total energy of primary particle smaller than the mass (e.g., an alpha particle with 1 MeV).

Best, Alvaro

andresailer commented 7 months ago

Doesn't this rather mean to check against insufficient total energy, and not check against the mass of the particle? For my education, what happens if one simulates the alpha with 1 MeV total energy?

I am also not sure if you refer to the particle.tbl file, which allows one to define additional particles, but it doesn't allow one to overwrite particles Geant4 already knows.

https://github.com/AIDASoft/DD4hep/blob/4db360e82ee16b5c822b7a3332e01ff225e5268f/DDG4/plugins/Geant4ExtraParticles.cpp#L110-L112

There is no way for us to know what the correct mass of a particle is. If there are MC records used as input, which define include the mass of particles, we could check that that mass is consistent with the mass of the defined, but this doesn't sound like the issue you are raising here.

Cheers, Andre

atolosadelgado commented 7 months ago

Hi,

For my education, what happens if one simulates the alpha with 1 MeV total energy?

ddsim still assigns a positive kinetic energy to the alpha, but well below 1 MeV. Simulation runs but the output is wrong. It can be very misleading.

we could check that that mass is consistent with the mass of the defined

Yes, that would be my point. Sorry I was not clear enough.

Best, Alvaro

andresailer commented 7 months ago

Yes, that would be my point. Sorry I was not clear enough.

But my proposal would do nothing for the 1 MeV Alpha case? The mass of the Alpha isn't changed, is it?

atolosadelgado commented 7 months ago

The alpha at 1 MeV should still trigger an exception because the tabulated mass is bigger than the input value of E_kin+mass*c**2.

andresailer commented 7 months ago

What is the tabulated mass? What is E_kin and mass in the 1 MeV alpha case? (I need to understand, or I can't even begin to implement anything to check this)

atolosadelgado commented 7 months ago

Geant4 maintains a table with all possible particles and their properties. Accessing the static properties of particles can be done as follows:

// Geant4 example: extended/electromagnetic/TestEm18/src/PrimaryGeneratorAction.cc
#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"
...
G4ParticleDefinition* particle = G4ParticleTable::GetParticleTable()->FindParticle("alpha");
G4double particle_mass = particle->GetPDGMass();
G4bool particle_stable   = particle->GetPDGStable();

The mass of the alpha particle is retrieved using the above snippet.

The DD4hep native particle gun requires the energy parameter to include the mass of the particle, so specifying 1 MeV for an alpha particle is incorrect. If the mass is not accounted for, it could lead to erroneous simulations.

Similar issues may arise when setting up the primaries from a hepmc3 file.

In my opinion, it would be safer to throw an exception when encountering such inconsistencies before starting the actual simulation.

There is a dedicated way to declare new particles for Geant4; it is not possible to run a simulation with an incorrect mass for a known particle.

MarkusFrankATcernch commented 5 months ago

Hi Alvaro,

I verified again in the code:

I cannot see anything wrong with that - except that the property name is perhaps not entirely detailed. What do I miss?

It is true, if an input file contains 3-momentum and energy and at the same time the particle mass loaded from external file is largely bigger than what was used to generate the input file, this might end up with a computed "negative" mass or momentum. Is this what you meant?

atolosadelgado commented 5 months ago

It is true, if an input file contains 3-momentum and energy and at the same time the particle mass loaded from external file is largely bigger than what was used to generate the input file, this might end up with a computed "negative" mass or momentum. Is this what you meant?

Yes, that it is. Additionally, even with positive mass, there might be a mismatch between the derived mass from input 3momentum+energy and the mass value used by Geant4. For example, if the input is E=5, pc=4, the derived mass would be mc2 = 3, but the internal Geant4 value for that particle mass may be 2.5.

In case the derived mass value is smaller than 10 eV, the dynamical mass is set to zero. In case is larger by more than 10 eV than the PDG value, the input mass is taken. This behavior is implemented here.

Thank you for taking a look to this :)

atolosadelgado commented 5 months ago

I had a quick chat with some Geant4 people, and I will open a review of the G4PrimaryParticle class to add a sanity check there.

A piece of code may be added in the hepMC3 reader in order to check the mass consistency of the derived value from the input and the value from Geant4/PDG. What do you think?

MarkusFrankATcernch commented 5 months ago

Ideally the check is done in a central piece of code where the data of all readers, generators and guns are processed. I suspect this is somewhere in Geant4InputHandling.cpp.

andresailer commented 5 months ago

Maybe here https://github.com/AIDASoft/DD4hep/blob/58fa8c0b16458b645a3edfd8469fddd87d8d8689/DDG4/src/Geant4InputHandling.cpp#L346

MarkusFrankATcernch commented 5 months ago

@andresailer This is a very good place. all is present and also the particle definitions are queried from Geant4. I will put the fix there.

atolosadelgado commented 5 months ago

after reading a bit more Geant4 code, I noticed that G4PrimaryParticle actually implements a check for negative masses here.

I have asked some colleagues from the FCC physics performance studies group, if implementing a check for the input particle mass would be an obstacle for their work or not (just to be sure).

MarkusFrankATcernch commented 5 months ago

Please see MR https://github.com/AIDASoft/DD4hep/pull/1259