Geant4PrimaryHandling: fix issue with multiple vertices in Geant4 GeneralParticleSource

DDG4/include/DDG4/Geant4InputHandling.h

@@ -42,7 +42,7 @@ namespace dd4hep {
     Geant4Particle* createPrimary(int particle_id, const Geant4Vertex* v, const G4PrimaryParticle* g4p);
 
     /// Create a DDG4 interaction record from a Geant4 interaction defined by a primary vertex
-    Geant4PrimaryInteraction* createPrimary(int mask, Geant4PrimaryMap* pm, const G4PrimaryVertex* gv);
+    Geant4PrimaryInteraction* createPrimary(int mask, Geant4PrimaryMap* pm, std::set<G4PrimaryVertex*>const& primaries);
 
     /// Initialize the generation of one event
     int generationInitialization(const Geant4Action* caller,const Geant4Context* context);

DDG4/src/Geant4GeneratorWrapper.cpp

@@ -70,18 +70,14 @@ void Geant4GeneratorWrapper::operator()(G4Event* event)  {
   Geant4PrimaryMap*   primaryMap = context()->event().extension<Geant4PrimaryMap>();
   set<G4PrimaryVertex*> primaries;
 
+  // Now generate the new interaction
+  generator()->GeneratePrimaryVertex(event);
+
   /// Collect all existing interactions (primary vertices)
   for(G4PrimaryVertex* v=event->GetPrimaryVertex(); v; v=v->GetNext())
     primaries.insert(v);
 
-  // Now generate the new interaction
-  generator()->GeneratePrimaryVertex(event);
-
-  /// Add all the missing interactions (primary vertices) to the primary event record.
-  for(G4PrimaryVertex* gv=event->GetPrimaryVertex(); gv; gv=gv->GetNext())  {
-    if ( primaries.find(gv) == primaries.end() )   {
-      Geant4PrimaryInteraction* inter = createPrimary(m_mask, primaryMap, gv);
-      prim->add(m_mask, inter);
-    }
-  }
+  // Add all the missing interactions (primary vertices) to the primary event record.
+  Geant4PrimaryInteraction* inter = createPrimary(m_mask, primaryMap, primaries);
+  prim->add(m_mask, inter);
 }

DDG4/src/Geant4InputHandling.cpp

@@ -88,6 +88,11 @@ static void collectPrimaries(Geant4PrimaryMap*         pm,
                              Geant4Vertex*             particle_origine,
                              G4PrimaryParticle*        gp)
 {
+  //if the particle is in the map, we do not have to do anything
+  if ( pm->get(gp) )  {
+    return;
+  }
+
   int pid = int(interaction->particles.size());
   Geant4Particle* p = createPrimary(pid,particle_origine,gp);
   G4PrimaryParticle* dau = gp->GetDaughter();
@@ -103,14 +108,13 @@ static void collectPrimaries(Geant4PrimaryMap*         pm,
 
   if ( dau )   {
     Geant4Vertex* dv = new Geant4Vertex(*particle_origine);
-    int vid = int(interaction->vertices.size());
     PropertyMask reason(p->reason);
     reason.set(G4PARTICLE_HAS_SECONDARIES);
 
     dv->mask = mask;
     dv->in.insert(p->id);
 
-    interaction->vertices[vid].emplace_back(dv) ;
+    interaction->vertices[mask].emplace_back(dv) ;
 
     for(; dau; dau = dau->GetNext())
       collectPrimaries(pm, interaction, dv, dau);
@@ -121,18 +125,19 @@ static void collectPrimaries(Geant4PrimaryMap*         pm,
 Geant4PrimaryInteraction* 
 dd4hep::sim::createPrimary(int mask,
                            Geant4PrimaryMap* pm,
-                           const G4PrimaryVertex* gv)
+                           std::set<G4PrimaryVertex*>const& primaries)
 {
   Geant4PrimaryInteraction* interaction = new Geant4PrimaryInteraction();
-  Geant4Vertex* v = createPrimary(gv);
-  int vid = int(interaction->vertices.size());
   interaction->locked = true;
   interaction->mask = mask;
-  v->mask = mask;
-  interaction->vertices[vid].emplace_back(v);
-
-  for (G4PrimaryParticle *gp = gv->GetPrimary(); gp; gp = gp->GetNext() )
-    collectPrimaries(pm, interaction, v, gp);
+  for (auto const& gv: primaries) {
+    Geant4Vertex* v = createPrimary(gv);
+    v->mask = mask;
+    interaction->vertices[mask].emplace_back(v);
+    for (G4PrimaryParticle *gp = gv->GetPrimary(); gp; gp = gp->GetNext()) {
+      collectPrimaries(pm, interaction, v, gp);
+    }
+  }
   return interaction;
 }
 
@@ -180,13 +185,14 @@ static void appendInteraction(const Geant4Action* caller,
   }
   Geant4PrimaryInteraction::VertexMap::iterator ivfnd, iv, ivend;
   for( iv=input->vertices.begin(), ivend=input->vertices.end(); iv != ivend; ++iv )   {
-    ivfnd = output->vertices.find((*iv).first) ; //(*iv).second->mask);
+    int theMask = input->mask;
+    ivfnd = output->vertices.find(theMask);
     if ( ivfnd != output->vertices.end() )   {
       caller->abortRun("Duplicate primary interaction identifier!",
                        "Cannot handle 2 interactions with identical identifiers!");
     }
     for(Geant4Vertex* vtx :  (*iv).second )
-      output->vertices[(*iv).first].emplace_back( vtx->addRef() );
+      output->vertices[theMask].emplace_back( vtx->addRef() );
   }
 }