Closed kkauder closed 1 year ago
@wdconinc
Hi @kkauder,
As far as I can tell it already goes wrong here
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?
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.
What about the event_pos()
you are using?
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.
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()?
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
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
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 :grin:
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
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;
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
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.
make install
cat < 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
ddsim --compactFile ./DD4hep/install/DDDetectors/compact/SiD.xml --numberOfEvents 1 --inputFiles small.hepmc --outputFile small.edm4hep.root
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 [...]
MCParticles.vertex.x = 0, 30, 30 MCParticles.vertex.y = 0, 0, 0 MCParticles.vertex.z = 0, 100, 100