Open simonge opened 8 months ago
Hi @simonge
Some kind of example (just the lines or a hepmc3 file) would be helpful.
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!
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://github.com/AIDASoft/DD4hep/blob/master/DDG4/hepmc/HepMC3EventReader.cpp
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.
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).
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.
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!
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.
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.
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?
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?
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...
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:
Any other suggestions or help would be appreciated.