Final state particles in HepMC are forced to vertex 0,0,0

[kkauder]

We are trying to originate particles from a specific vertex (for background studies). Since HepMC needs an incoming particle in every vertex, we set this minimal test up like so:

g (13) g(1)
----> v --->

All particles are gammas, the number in parentheses is their status. 13 is arbitrary to signify "virtual", anything above 10 is ok for HepMC. v is at (30,0,100,100). In the example below I add a second out-going gamma with status 13 to demonstrate that it has the correct vertex.

• OS version: Some debian flavor (I'm running in a container)
• Compiler version: GCC 12.2.0
• Package version: latest commit in git log is 665d9de
• Reproduced by:

cmake -S DD4hep -B DD4hep/build -DDD4HEP_USE_GEANT4=ON -DDD4HEP_USE_EDM4HEP=ON -DDD4HEP_USE_HEPMC3=ON -DBoost_NO_BOOST_CMAKE=ON -DDD4HEP_USE_LCIO=ON -DBUILD_TESTING=ON -DROOT_DIR=/usr/local -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=DD4hep/install

make install

Make an input file

cat <<EOF > small.hepmc
HepMC::Version 3.02.00
HepMC::Asciiv3-START_EVENT_LISTING
E 1 1 3
U GEV MM
P 1 0 22 1.4367394278317271e+00 0.0000000000000000e+00 4.7891314261057571e+00 5.0000000000000000e+00 0.0000000000000000e+00 13
V -1 0 [1] @ 3.0000000000000000e+01 0.0000000000000000e+00 1.0000000000000000e+02 1.0000000000000000e+02
P 2 -1 22 5.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 5.0000000000000000e+00 0.0000000000000000e+00 1
P 3 -1 22 10.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 10.0000000000000000e+00 0.0000000000000000e+00 13
EOF

Then run with

ddsim --compactFile ./DD4hep/install/DDDetectors/compact/SiD.xml --numberOfEvents 1 --inputFiles small.hepmc --outputFile small.edm4hep.root

Result:

root -l small.edm4hep.root
events->Show(0)
[...]
 MCParticles.PDG = 22, 22, 22
 MCParticles.generatorStatus = 13, 1, 13
[...]
 MCParticles.vertex.x = 0, 0, 30
 MCParticles.vertex.y = 0, 0, 0
 MCParticles.vertex.z = 0, 0, 100
[...]

• Input: see above. Quicker to just copy paste
• Output: I'll be happy to update if you really want the full thing. There's a lot of cmake and geant happening...
• Goal: A short description of what you are trying to achieve
  Both out-going particles, but especially the final one that I actually care about, should start from the vertex. I.e.

 MCParticles.vertex.x = 0, 30, 30
 MCParticles.vertex.y = 0, 0, 0
 MCParticles.vertex.z = 0, 100, 100

Notes:

• This is the shortest demo. I've also constructed more complete examples with proper beam (status 4) particles to no effect.
• I've chased it down further into Geant4InputHandling.cpp, Geant4InputAction.cpp, it gets increasingly difficult for me to figure out what's wrong.
• I've done some work in the tedium that is reading HepMC in a different package, including displaced vertices and pre-assigned decays
  https://github.com/eic/east/blob/main/PrimGenInterface/src/eASTHepMC3Interface.cc
  I'll be happy to help if you point me in the right direction

[kkauder]

@wdconinc

[andresailer]

Hi @kkauder,

As far as I can tell it already goes wrong here

DD4hep/DDG4/hepmc/HepMC3EventReader.cpp

Lines 132 to 142 in cc4c4c5

	if ( p->parents.size() == 0 ) {
	Geant4Vertex* vtx = new Geant4Vertex ;
	vertices.emplace_back( vtx );
	vtx->x = p->vsx;
	vtx->y = p->vsy;
	vtx->z = p->vsz;
	vtx->time = p->time;
	vtx->out.insert(p->id) ;
	}

I tried printing also genEvent.event_pos()
like you are using here
https://github.com/eic/east/blob/7821fb87b337de5ecb7f526ef8b45102b5f3da73/PrimGenInterface/src/eASTHepMC3Interface.cc#L128-L129

  auto pos = hepevt.event_pos();
  auto* g4vtx  = new G4PrimaryVertex(pos.x()*mm, pos.y()*mm, pos.z()*mm, 0);

But that also gives me 0 0 0. There is nothing in dd4hep that manipulates the genEvent, so the information is not (yet) in the hepmc::genevent?
Is there some step in your code I missed that moves the vertex position?

[kkauder]

That part shouldn't be the one that matters. It only kicks in for parent-less particles, i.e. the first virtual gamma. I care about the one that does have a parent, and thus a proper production vertex.

[andresailer]

What about the event_pos() you are using?

[kkauder]

Good question, I don't think event_pos() is quite appropriate for this fragment of an event. I'd probably define that as the vertex where the beam particles meet.
But my point is, the gamma I care about has a clearly defined origin vertex (and its virtual sister even remembers it), that gets lost somewhere.

[andresailer]

What does HepMC define as the event_pos though? and where do you set the vertex position in your eAST hepmc3 interface if not via event_pos()?

[kkauder]

I have to admit I don't know. The way I described it in the comments, event_pos() should have picked up the only vertex in the event, not 0, 0, 0

[kkauder]

From https://gitlab.cern.ch/hepmc/HepMC3/-/blob/master/src/GenEvent.cc

const FourVector& GenEvent::event_pos() const {
    return m_rootvertex->data().position;
}

Default ctor sets it to 0. Lookin g further, it seems to be a strange beast that is not written out. Can be manipulated with shift_position_by though.

That said, assuming I'd create the event by using this root vertex, that would mean
a) Only one shifted vertex is allowed (our events have many)
b) Preassigned decays still have to be treated differently since they don't come from the root vertex. So the origin vertex of final state particles has to be respected regardless of rootvertex

[andresailer]

1. Do you mean multiple primary vertices, or non-primary vertices from preassigned decays?
   Can one do multiple primary interactions in hepmc3?

2. In the DD4hep logic preassigned decays are treated differently: but all particles between the primary vertex and the preassigned decay must either be known to Geant4, or have zero-lifetime. Geant has to do the propagation to account for magnetic fields or material effects.
   This is why only the parents.size()==0 particles are used to created G4Vertices.

3. We will discuss on Monday, probably easier to go into more details

4. Might be easier to write an EDM4hep reader 😁

5. If the first vertex is always the interaction location, we could maybe just take the "end" position of the parentless particle as the G4Vertex position, but I have not enough understanding of the hepmc format or API

[kkauder]

After discussing in person with @andresailer, his proposed solution works!

# git diff DDG4/hepmc/HepMC3EventReader.cpp
diff --git a/DDG4/hepmc/HepMC3EventReader.cpp b/DDG4/hepmc/HepMC3EventReader.cpp
index 67c576ff..2236ce69 100644
--- a/DDG4/hepmc/HepMC3EventReader.cpp
+++ b/DDG4/hepmc/HepMC3EventReader.cpp
@@ -133,9 +133,9 @@ HEPMC3EventReader::readParticles(int event_number, Vertices& vertices, Particles

       Geant4Vertex* vtx = new Geant4Vertex ;
       vertices.emplace_back( vtx );
-      vtx->x = p->vsx;
-      vtx->y = p->vsy;
-      vtx->z = p->vsz;
+      vtx->x = p->vex;
+      vtx->y = p->vey;
+      vtx->z = p->vez;
       vtx->time = p->time;

[kkauder]

I'm working on putting this into a PR but there's a snag. If a particle has no parent and no end vertex (such as a particle gun),

HepMC::Version 3.02.00
HepMC::Asciiv3-START_EVENT_LISTING
E 1 0 1
U GEV MM
P 1 -1 22 5.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 5.0000000000000000e+00 0.0000000000000000e+00 1

the code doesn't break but also doesn't do the right thing, instead falling back to 0,0,0 again.
In fact, it cannot do anything else since the production vertex isn't part of the P line

[kkauder]

#1017