Persist hepmc vertex status code into output.

[simonge]

Hepmc3 vertices contain a status code which there currently appears no way to keep to the DD4hep output (edm4hep).

We are wanting to use hepmc events which will contain several vertices coming from different processes/generators. To keep track of the process responsible for a vertex the hepmc3 status code can be used, but once the simulation is run this information is lost.

Potential solutions we have considered:

• Make use of some unused simulatorStatus bits in MCParticles.
• Save a (currently non-existent) MCVertex data type based on the hepmc3 vertex. In addition, if incoming particles aren't relevant or available the MCVertex would provide a reference point to associate the particles back to the process and the other associated particles.

Any other suggestions or help would be appreciated.

[andresailer]

Hi @simonge

Some kind of example (just the lines or a hepmc3 file) would be helpful.

[simonge]

Hi @andresailer,

This is a relatively simple example where a low-Q2 physics reaction is given the status code 1 and Bremsstrahlung vertices have status code 2. (I have stripped it down a bit for clarity)

E 0 23 93
U GEV MM
P 1 0 11 [px,py,px,E,m] 4
P 2 0 2212 [px,py,px,E,m] 4
V -1 **1** [1,2] @ [x,y,z,t]
P 3 -1 11 [px,py,px,E,m] 4
P 4 -1 2212 [px,py,px,E,m] 4
P 5 -1 11 [px,py,px,E,m] 1
P 6 0 11 [px,py,px,E,m] 4
P 7 0 2212 [px,py,px,E,m] 4
V -2 **2** [6,7] @ [x,y,z,t]
P 8 -2 22 [px,py,px,E,m] 1
P 9 -2 11 [px,py,px,E,m] 1
...
P 90 0 11  [px,py,px,E,m] 4
P 91 0 2212  [px,py,px,E,m] 4
V -23 **2** [90,91] @ [x,y,z,t]
P 92 -23 22  [px,py,px,E,m] 1
P 93 -23 11 [px,py,px,E,m] 1
E 2 22 89
U GEV MM
P 1 0 11 [px,py,px,E,m] 4
P 2 0 2212 [px,py,px,E,m] 4
V -1 **1** [1,2] @ [x,y,z,t]
P 3 -1 11 [px,py,px,E,m] 4
P 4 -1 2212 [px,py,px,E,m] 4
P 5 -1 11 [px,py,px,E,m] 1
P 6 0 11 [px,py,px,E,m] 4
P 7 0 2212 [px,py,px,E,m] 4
V -2 **2** [6,7] @ [x,y,z,t]
P 8 -2 22 [px,py,px,E,m] 1
P 9 -2 11 [px,py,px,E,m] 1
...

In Geant4EventReaderHepMC.cpp L502 the "status code" is read in as dummy.

Thanks!

[andresailer]

HI @simonge

Thanks!

Don't look at the Geant4EventReaderHepMC, that is not fully working. You really want to build against HepMC3 and look at https:/AIDASoft/DD4hep/blob/master/DDG4/hepmc/HepMC3EventReader.cpp

[wdconinc]

This was brought up in today's key4hep meeting, and the main limitation here is that the MCParticle data type in EDM4hep does not have any field for storing the vertex ID. As far as I can tell, it is the only field in the HepMC3 input event records that cannot be mapped onto MCParticle cleanly.

Reading mcp->production_vertex()->status() gets the status into DD4hep, but there is not really a place to write it out.

[tmadlener]

Could some of the free bits in the generatorStatus of the edm4hep::MCParticle be used for this? Obviously slightly contingent on how many bits would actually be needed to fully map this.

Edit: Looking at the edm4hep.yaml file it looks like there is no reserved bits yet in generatorStatus (TBC).

[wdconinc]

That was our idea: to reuse some unused bits. Nothing immediately obvious jumps out.

The alternative (which would solve some other requests we've had over the past years, and was also included in the original issue by Simon) is to add a MCVertex data type in addition to MCParticle. That would be light-weight (just position and vertex status are included HepMC3, and then relations to in and out particles), and could be used to store step-wise trajectories as well if we extend it to beyond the hard scattering graph. That approach would be backwards compatible since it wouldn't require repurposing bits.

[MarkusFrankATcernch]

These are the bits avaible in the MCParticle:
MCParticle::status Bit 18 -> 31 (including) of the status word can be used by users.
MCParticle::genStatus Bit 5 to 15 (including) of the status word can be used by users.
MCParticle::_spare Bit 0 to Bit 7 (including): word may be repurposed.

I think the MCParticle::genStatus is probably the most appropriate.
Important:
If you use the bits for HepMC3 please update the documentation in the header!

[simonge]

I was hesitant to suggest MCParticle::genStatus as it has the potential to break and complicate a lot of user code. Where at the moment MCParticle::genStatus==1 means the same thing in HepMC3 and the edm4hep output, mixing in the vertex status code will make this much less clear.

[MarkusFrankATcernch]

If statements like MCParticle::genStatus==1 are used this is a mistake. These fields were always meant to be bifields....
We should try to identify them and fix these.

[andresailer]

One issue for genStatus is that people might select in ROOT tree based analysis on that number, where using a bitfield encoder is non-trivial?

[wdconinc]

Would it be relatively straightforward to implement this as an option (default off), or is that suboptimal because of the additional config and code paths required and a different class of errors?

[gaede]

To add my 2 cents: in LCIO and therefore also in EDM4hep we are using the MCParticle::genStatus as int field and not as bit field, ie. mcp.getGeneratorStatus()==1 is perfectly valid code and used in very many places. So whatever is done with the dd4hep::MCParticle::genStat bits - when written to LCIO and EDM4hep etxtra bits should be scrapped. One could use bits in MCParticle::simStatus as this is a true bitfield. Another option would be to add an MCParticleVertex object to EDM4hep - maybe a reasonable extension that would also make conversions from HepMC easier. Would need some discussions in EDM4hep though...