Vertex time not imported correctly from hepmc

[kkauder]

!!!IMPORTANT!!!: To facilitate faster and easier response to your issue please provide in addition to the description of the issue also the following information

• OS version: e.g. centos7
  container based on debian,
  https://github.com/eic/eic-shell
  https://github.com/eic/containers

• Compiler version: e.g GCC 8.0
  gcc 12

• Package version: tag or commit ID, or GitHub branch
  dd4hep@1.24 +ddg4 +ddcad -frames +hepmc3 +lcio
  (None of the above matters. We are stuck for technical reasons with 1.24 now, but the current branch has the same code.)

• Reproduced by: exact steps to reproduce the problem: checkout, setup environment (Geant4/ROOT version, LCG/iLCSoft/Key4hep release), cmake, build, run...
  Run the attached file through any detector

ddsim -N 1 --compactFile $DETECTOR_PATH/$DETECTOR_CONFIG.xml --inputFile simple.hepmc --outputFile simple.edm4hep.root

It has three final photons emerging from different vertices and crucially at different times, 0, 1000, and 2000 ns

• Input:
  simple.hepmc.zip

HepMC::Version 3.02.05
HepMC::Asciiv3-START_EVENT_LISTING
E 1 3 6
U GEV MM
P 1 0 22 5.4900591458096715e-08 -2.3360014308601191e-09 -6.3873918346503390e-06 6.3876281968507648e-06 0.0000000000000000e+00 4
V -1 0 [1] @ 2.9390599999999999e+01 -1.2505599999999999e+00 -3.4194400000000001e+03 0.0000000000000000e+05
P 2 -1 22 2.5091952042300000e-08 -4.2977507671200002e-10 -6.3875788989600000e-06 6.3876281968507648e-06 0.0000000000000000e+00 1
P 3 0 22 -2.2529055020399048e-07 1.7158543817213154e-08 -9.4934358251372510e-06 9.4961241574463718e-06 0.0000000000000000e+00 4
V -2 0 [3] @ -2.9304700000000000e+01 2.2319000000000000e+00 -1.2348599999999999e+03 3.0000000000000000e+05
P 4 -2 22 -6.0711923840800009e-08 1.9305896843600003e-10 -9.4959300776000006e-06 9.4961241574463718e-06 0.0000000000000000e+00 1
P 5 0 22 3.5707443339668777e-08 3.1812558302218797e-09 -5.3356106857103988e-06 5.3357311149799803e-06 0.0000000000000000e+00 4
V -3 0 [5] @ 2.9271799999999999e+01 2.6078899999999998e+00 -4.3739600000000000e+03 6.0000000000000000e+05
P 6 -3 22 2.1607412136099995e-08 -7.3254237170000003e-10 -5.3356873141600002e-06 5.3357311149799803e-06 0.0000000000000000e+00 1
HepMC::Asciiv3-END_EVENT_LISTING

• Output: full build and run log output
  You don't want that... Let me zoom in on the problem, close to the end

PrimaryHandler   INFO  +++++ G4PrimaryVertex at (+2.94e+01,-1.25e+00,-3.42e+03) [mm] +0.00e+00 [ns]
PrimaryHandler   INFO  +++++ G4PrimaryVertex at (-2.93e+01,+2.23e+00,-1.23e+03) [mm] +0.00e+00 [ns]
PrimaryHandler   INFO  +++++ G4PrimaryVertex at (+2.93e+01,+2.61e+00,-4.37e+03) [mm] +0.00e+00 [ns]

• Goal: A short description of what you are trying to achieve

PrimaryHandler   INFO  +++++ G4PrimaryVertex at (+2.94e+01,-1.25e+00,-3.42e+03) [mm] +0.00e+00 [ns]
PrimaryHandler   INFO  +++++ G4PrimaryVertex at (-2.93e+01,+2.23e+00,-1.23e+03) [mm] +1.00e+03 [ns]
PrimaryHandler   INFO  +++++ G4PrimaryVertex at (+2.93e+01,+2.61e+00,-4.37e+03) [mm] +2.00e+03 [ns]

The problem is here

DD4hep/DDG4/hepmc/HepMC3EventReader.cpp

Lines 104 to 105 in 504cc12

	p->time = vsx.get_component(3) * len_unit / CLHEP::c_light;
	p->properTime = vsx.get_component(3) * len_unit / CLHEP::c_light;

and here:

DD4hep/DDG4/hepmc/HepMC3EventReader.cpp

Lines 136 to 149 in 504cc12

	if ( p->parents.size() == 0 ) {
	// A particle without a parent in HepMC3 can only be (something like) a beam particle, and it is attached to the
	// root vertex, by default (0,0,0) and equal for all parent-less particles. Therefore we can take the end vertex
	// of the parentless particle as the start vertex for outgoing particles. Note that for a particle without end
	// vertex (such as in a particle gun), it defaults to (0,0,0). This cannot be fixed, the information simply isn't
	// in the HepMC file. Having a parent enforces a vertex, having no parent forbids a vertex.
	Geant4Vertex* vtx = new Geant4Vertex ;
	vertices.emplace_back( vtx );
	vtx->x = p->vex;
	vtx->y = p->vey;
	vtx->z = p->vez;
	vtx->time = p->time;

Particles correctly know their start vertex and, from there, creation time. However, GEANT vertices are only constructed from the end vertices of particles without a parent. I have a simple patch for that at https://github.com/eic/eic-spack/blob/develop/packages/dd4hep/vertex-time.patch

The reason I'm not making a pull request from that is that it may be better to rethink the logic of using parent-less particles for geant vertices. I have some idea about why it's done but I don't think it's right. And I believe it should also not work for preassigned decays. I have some experience with reading HepMC into GEANT, https://github.com/eic/east/blob/main/PrimGenInterface/src/eASTHepMC3Interface.cc and I'd be happy to help. Otherwise please let me know and I'll open a PR for the kludge.

[andresailer]

Why don't you make a PR with the patch to fix this problem?

[kkauder]

It's hacky. Because it adds a data member vet (vertex end time) to Geant4Particle. In principle that should also propagate through the getters/setters and I haven' done that. It simply compounds the current weird logic with some glue. But I will do it.

[kkauder]

Keep in mind, this won't fix secondaries / pre-assigned decays

[andresailer]

Can you elaborate on the problem of secondaries and pre assigned decays there is?

[kkauder]

The code as it is only generates GEANT vertices for parentless particles, which e.g. a decayed J/psi and its daughters cannot be.

[andresailer]

Do you have a hepmc file that shows that there actually is a problem? Because pre-assigned decays are working right now, for example in this file https://github.com/AIDASoft/DD4hep/blob/master/DDTest/inputFiles/Pythia_output.hepmc

[andresailer]

run ddsim --output.inputStage VERBOSE [ETC]
and you see

Vertex  ( 1.41881e-12[mm], -1.93334e-12[mm], 0[mm], 0[ns] ) Weight 1
  -- Primary particles ::    # of primaries =37
==== PDGcode -211  Particle name pi-
 Assigned charge : -1
     Momentum ( 0.218022[GeV/c], -0.209761[GeV/c], 8.18433[GeV/c] )
     kinetic Energy : 8.05154 [GeV]
     Mass : 0.13957 [GeV]
     Polarization ( 0, 0, 0 )
     Weight : 1
==== PDGcode 111  Particle name pi0
 Assigned charge : 0
     Momentum ( -0.445784[GeV/c], -0.865636[GeV/c], 9.7899[GeV/c] )
     kinetic Energy : 9.70415 [GeV]
     Mass : 0.13498 [GeV]
     Polarization ( 0, 0, 0 )
     Weight : 1
     PreAssigned proper decay time : 2.27179e-08 [ns] 
>>>> Daughters
==== PDGcode 22  Particle name gamma
 Assigned charge : 0
     Momentum ( -0.0192403[GeV/c], -0.00156629[GeV/c], 0.209193[GeV/c] )
     kinetic Energy : 0.210081 [GeV]
     Mass : 0 [GeV]
     Polarization ( 0, 0, 0 )
     Weight : 1
==== PDGcode 22  Particle name gamma
 Assigned charge : 0
     Momentum ( -0.426543[GeV/c], -0.86407[GeV/c], 9.58071[GeV/c] )
     kinetic Energy : 9.62905 [GeV]
     Mass : 0 [GeV]
     Polarization ( 0, 0, 0 )
     Weight : 1
<<<< End of link
...

[kkauder]

Interesting, let me give that a try. I made assumptions I probably shouldn't have.

[andresailer]

Pre-assigned decays do not need new G4Vertices, they need the decay products to be assigned to particles:

Which is done here:

DD4hep/DDG4/src/Geant4InputHandling.cpp

Line 430 in 8c5fd05

p4->SetDaughter((*i).second);

[kkauder]

You are correct, I withdraw my comment :) I'm currently getting

RootOutput       ERROR +++ [Event:0] Exception while saving event:std::bad_alloc
RootOutput       ERROR +++ [Event:0] Exception while saving event:std::bad_alloc
Geant4Kernel     FATAL +++ Exception while simulating:std::bad_alloc

when running over this pythia event, but I suspect that's because it's going through the full detector and it exceeds the RAM allotted to my docker image somewhere in the writing stage. I'll follow up if it persists after choosing a simpler geometry.

Your PR works, and after it's merged I consider this issue fully resolved. Thanks a lot for the fast and smart help!!

[andresailer]

Yes, that is a very heavy event. I usually disable almost everything in a detector when debugging MC record issues, and it still takes minutes to simulate something from this file.

Cheers!

[kkauder]

To confirm, if I use just the inner tracker it works just fine. Thanks again!