ddsim does not simulate when pythia vertex smearing is enabled

[scott-snyder]

hi -

I've been running into issues with ddsim crashing with errors like

Input            FATAL +++ Stable particle (PDG: -16       ) with daughters! check your MC record, adapt particle.tbl file...

when run with vertex smearing.

This seems to be fallout from this commit:

commit eb00b5bf7035f6867f0369f25dc187f2852002f8
Author: Andre Sailer <andre.philippe.sailer@cern.ch>
Date:   Tue Jul 26 10:51:24 2022 +0200

    Geant4InputHandling: add sanity check for stable particle with children, end simulation

As an example, this was executed on a machine running

LSB Version:    :core-4.1-amd64:core-4.1-ia32:core-4.1-noarch
Distributor ID: Scientific
Description:    Scientific Linux release 7.3 (Nitrogen)
Release:        7.3
Codename:       Nitrogen

with the current key4hep setup (which gives gcc 11.2.0 and DD4hep 1.23).

Here is an example reproducer. The whizard input files are attached
(though you have to rename them from .txt because githib chokes on them otherwise).
They generate 500 events of (Z->mumu)(H->invisible), with whizard's
built-in pythia enabled, configured with the tune used for FCC generation.
In the first file, hinvmm1.sin, vertex smearing is enabled in pythia
as in the FCC tune. The second, hinvmm2.sin, is the same except that the
vertex smearing has been disabled.

$ . /cvmfs/sw.hsf.org/key4hep/setup.sh
# Currently gives key4hep-stack/2023-04-08
$ mkdir hinvmm1 hinvmm2
$ cd hinvmm1
$ whizard -L hinvmm1.log ../hinvmm1.sin 2>&1|tee log
$ cd ../hinvmm2
$ whizard -L hinvmm2.log ../hinvmm2.sin 2>&1|tee log
$ cd ..
$ git clone https://github.com/iLCSoft/CLICPerformance
$ ddsim --compactFile $LCGEO/FCCee/compact/FCCee_o1_v05/FCCee_o1_v05.xml --inputFiles hinvmm1/hinvmm1.stdhep --outputFile hinvmm1.sim.edm4hep.root --numberOfEvents=100  --steeringFile CLICPerformance/fcceeConfig/fcc_steer.py 2>&1|tee log-hinvmm1

This fails with

PrimaryHandler   INFO  +++++ G4PrimaryVertex at (+0.00e+00,+0.00e+00,+0.00e+00) [mm] +0.00e+00 [ns]
Input            FATAL +++ Stable particle (PDG: -16       ) with daughters! check your MC record, adapt particle.tbl file...
Geant4Kernel     FATAL +++ Exception while simulating:Cannot Simmulate this MC Record

I also tried using ddsim's smearing directly rather than having
pythia do it. This does not work either:

$ ddsim --compactFile $LCGEO/FCCee/compact/FCCee_o1_v05/FCCee_o1_v05.xml --inputFiles hinvmm2/hinvmm2.stdhep --outputFile hinvmm2.sim.edm4hep.root --numberOfEvents=100   --vertexSigma 0.0098 2.54e-5 0.646 1.9372  --steeringFile CLICPerformance/fcceeConfig/fcc_steer.py 2>&1|tee log-hinvmm2

Fails with

Geant4InteractionVertexSmear INFO  +++ Smearing primary vertex for interaction type 4 (1 Vertices, 21 particles) by (-6.03e-03 mm, +2.29e-05 mm, +1.93e-01 mm, +2.62e+00 ns)
PrimaryHandler   INFO  +++++ G4PrimaryVertex at (-6.03e-03,+2.29e-05,+1.93e-01) [mm] +2.62e+00 [ns]
Input            FATAL +++ Stable particle (PDG: 12        ) with daughters! check your MC record, adapt particle.tbl file...
Geant4Kernel     FATAL +++ Exception while simulating:Cannot Simmulate this MC Record

The issue seems to stem at least partly from numerical issues where
neutrino masses calculated as sqrt(E^2 - p^2) are not exactly zero.
I found that clipping very small masses to zero solves the issues
i was seeing, but this doesn't seem like the best way to fix this.

diff --git a/DDG4/src/Geant4InputHandling.cpp b/DDG4/src/Geant4InputHandling.cpp
index 88248710..c9e9066e 100644
--- a/DDG4/src/Geant4InputHandling.cpp
+++ b/DDG4/src/Geant4InputHandling.cpp
@@ -380,7 +380,8 @@ getRelevant(set<int>& visited,
     int first_daughter = *(dau.begin());
     Geant4ParticleHandle dp = pm[first_daughter];
     double en = p.energy();
-    double me = en > std::numeric_limits<double>::epsilon() ? p->mass / en : 0.0;
+    double ma = p->mass < 1e-2 ? 0 : p->mass;
+    double me = en > std::numeric_limits<double>::epsilon() ? ma / en : 0.0;
     //  fix by S.Morozov for real != 0
     double proper_time = fabs(dp->time-p->time) * me;
     double proper_time_Precision = pow(10.,-DBL_DIG)*me*fmax(fabs(p->time),fabs(dp->time));

hinvmm1.sin.txt
hinvmm2.sin.txt

[andresailer]

ddsim doesn't crash, it refuses to simulate. I think the code works as intended. (Without the check it would probably result in infinite geant4 simulation or an actual crash... It just runs, but doesn't put the neutrino in the right starting point)

The problem is that there actually is a stable particle with a daughter...
(I ran stdhepjob to convert the stdhep file to slcio to be able to dump it directly, lcio output for whizard also works.)

[00000016]   11|        14| 1.91e+01,-5.26e+00, 1.60e+01| 0.00e+00, 0.00e+00, 0.00e+00| 2.55e+01| 2 |[    0   ]| 0.00e+00, 0.00e+00, 0.00e+00|-1.13e-02, 1.56e-05, 1.28e+00| 6.74e-07| 0.00e+00| 0.00e+00, 0.00e+00, 0.00e+00|  (0, 0)   | [6] - [17] 
[snip]
[00000022]   17|        14| 1.91e+01,-5.26e+00, 1.60e+01| 0.00e+00, 0.00e+00, 0.00e+00| 2.55e+01| 1 |[    0   ]|-1.13e-02, 1.56e-05, 1.28e+00| 0.00e+00, 0.00e+00, 0.00e+00| 4.77e-07| 0.00e+00| 0.00e+00, 0.00e+00, 0.00e+00|  (0, 0)   | [11] - [] 

And according to the information in LCIO the first neutrino with status 2 flies a non-zero distance, which corresponds to a certain lifetime. So it is not the mass that is the problem, this particle is claimed to travel and then "decay".

Can you use the DD4hep vertex smear plugin instead, via the ddsim steering file, or via command line.

## FourVector of translation for the Smearing of the Vertex position: x y z t
SIM.vertexOffset = [0.0, 0.0, 0.0, 0.0]
## FourVector of the Sigma for the Smearing of the Vertex position: x y z t
SIM.vertexSigma = [0.0, 0.0, 0.0, 0.0]

[scott-snyder]

hi -

The second example i give uses the vertex smearing provided by ddsim; it fails in a similar manner.

[andresailer]

In the second example, why do you forward the event to pythia? it doesn't do anything but repeat the entries, does it?

[andresailer]

With #1097, the second example works for me, even when the particles are "duplicate".
The problem was the negative mass.

[scott-snyder]

Thanks! I'll give that a try.

Yes, using pythia is superfluous for Z->mumu, but we're also interested in Z->qq, and i'm trying to have a consistent setup, based on configurations in https://github.com/HEP-FCC/FCC-config, but changed to have whizard force an H->invisible decay. I sent Z->mumu as an example since that's the simplest and runs fastest.

[scott-snyder]

I can confirm that your change fixes this problem for all the final states we're looking at. Thanks!