Enable running ExternalDecays (photospp only for the moment) in the HepMC3 mode

GeneratorInterface/Core/interface/GeneratorFilter.h

@@ -183,10 +183,18 @@ namespace edm {
       // The external decay driver is being added to the system,
       // it should be called here
       //
-      if (decayer_) {  // handle only HepMC2 for the moment
-        auto t = decayer_->decay(event.get());
-        if (t != event.get()) {
-          event.reset(t);
+      if (decayer_) {
+        if (ivhepmc == 2) {  // handle HepMC2
+          auto t = decayer_->decay(event.get());
+          if (t != event.get()) {
+            event.reset(t);
+          }
+        }
+        if (ivhepmc == 3) {  // handle HepMC3
+          auto t = decayer_->decay(event3.get());
+          if (t != event3.get()) {
+            event3.reset(t);
+          }
         }
       }
       if (ivhepmc == 2 && !event.get())

GeneratorInterface/Core/test/FailingGeneratorFilter.cc

@@ -119,6 +119,7 @@ namespace test {
     std::vector<std::string> sharedResources() const { return {}; }
 
     HepMC::GenEvent* decay(HepMC::GenEvent const*) { return nullptr; }
+    HepMC3::GenEvent* decay(HepMC3::GenEvent const*) { return nullptr; }
     void statistics() const {}
 
     void init(const edm::EventSetup&) const {}

GeneratorInterface/ExternalDecays/BuildFile.xml

@@ -4,6 +4,7 @@
 <use name="heppdt"/>
 <use name="clhep"/>
 <use name="hepmc"/>
+<use name="hepmc3"/>
 <use name="GeneratorInterface/Core"/>
 <use name="GeneratorInterface/LHEInterface"/>
 <use name="GeneratorInterface/EvtGenInterface"/>

GeneratorInterface/ExternalDecays/interface/ExternalDecayDriver.h

@@ -12,6 +12,10 @@ namespace HepMC {
   class GenEvent;
 }
 
+namespace HepMC3 {
+  class GenEvent;
+}
+
 namespace CLHEP {
   class HepRandomEngine;
 }
@@ -40,6 +44,9 @@ namespace gen {
     HepMC::GenEvent* decay(HepMC::GenEvent* evt);
     HepMC::GenEvent* decay(HepMC::GenEvent* evt, lhef::LHEEvent* lheEvent);
 
+    HepMC3::GenEvent* decay(HepMC3::GenEvent* evt);
+    HepMC3::GenEvent* decay(HepMC3::GenEvent* evt, lhef::LHEEvent* lheEvent);
+
     void statistics() const;
 
     void setRandomEngine(CLHEP::HepRandomEngine*);

GeneratorInterface/ExternalDecays/src/ExternalDecayDriver.cc

@@ -8,6 +8,7 @@
 #include "GeneratorInterface/PhotosInterface/interface/PhotosFactory.h"
 #include "GeneratorInterface/PhotosInterface/interface/PhotosInterfaceBase.h"
 #include "HepMC/GenEvent.h"
+#include "HepMC3/GenEvent.h"
 #include "FWCore/Concurrency/interface/SharedResourceNames.h"
 // LHE Run
 #include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h"
@@ -63,6 +64,12 @@ HepMC::GenEvent* ExternalDecayDriver::decay(HepMC::GenEvent* evt, lhef::LHEEvent
   return decay(evt);
 }
 
+HepMC3::GenEvent* ExternalDecayDriver::decay(HepMC3::GenEvent* evt, lhef::LHEEvent* lheEvent) {
+  if (fTauolaInterface)
+    fTauolaInterface->SetLHE(lheEvent);
+  return decay(evt);
+}
+
 HepMC::GenEvent* ExternalDecayDriver::decay(HepMC::GenEvent* evt) {
   if (!fIsInitialized)
     return evt;
@@ -88,6 +95,33 @@ HepMC::GenEvent* ExternalDecayDriver::decay(HepMC::GenEvent* evt) {
   return evt;
 }
 
+HepMC3::GenEvent* ExternalDecayDriver::decay(HepMC3::GenEvent* evt) {
+  if (!fIsInitialized)
+    return evt;
+
+#if 0
+  if (fEvtGenInterface) {
+    evt = fEvtGenInterface->decay(evt);
+    if (!evt)
+      return nullptr;
+  }
+
+  if (fTauolaInterface) {
+    evt = fTauolaInterface->decay(evt);
+    if (!evt)
+      return nullptr;
+  }
+#endif
+
+  if (fPhotosInterface) {
+    evt = fPhotosInterface->apply(evt);
+    if (!evt)
+      return nullptr;
+  }
+
+  return evt;
+}
+
 void ExternalDecayDriver::init(const edm::EventSetup& es) {
   if (fIsInitialized)
     return;

GeneratorInterface/PhotosInterface/BuildFile.xml

@@ -2,6 +2,7 @@
 <use name="FWCore/Framework"/>
 <use name="FWCore/PluginManager"/>
 <use name="hepmc"/>
+<use name="hepmc3"/>
 <use name="clhep"/>
 <export>
   <lib name="1"/>

GeneratorInterface/PhotosInterface/interface/PhotosInterfaceBase.h

@@ -4,6 +4,7 @@
 #include "FWCore/ParameterSet/interface/ParameterSet.h"
 #include "FWCore/Framework/interface/EventSetup.h"
 #include "HepMC/GenEvent.h"
+#include "HepMC3/GenEvent.h"
 #include <vector>
 #include "CLHEP/Random/RandomEngine.h"
 #include "HepMC/SimpleVector.h"
@@ -21,6 +22,7 @@ namespace gen {
     virtual void init() = 0;
     virtual const std::vector<std::string>& specialSettings() { return fSpecialSettings; }
     virtual HepMC::GenEvent* apply(HepMC::GenEvent* evt) { return evt; }
+    virtual HepMC3::GenEvent* apply(HepMC3::GenEvent* evt) { return evt; }
     virtual void avoidTauLeptonicDecays() = 0;
     ;
     virtual void configureOnlyFor(int) = 0;

GeneratorInterface/PhotosInterface/interface/PhotosppInterface.h

@@ -15,6 +15,10 @@ namespace HepMC {
   class GenVertex;
 }  // namespace HepMC
 
+namespace HepMC3 {
+  class GenEvent;
+}  // namespace HepMC3
+
 namespace gen {
   class PhotosppInterface : public PhotosInterfaceBase {
   public:
@@ -25,6 +29,7 @@ namespace gen {
     void init() override;
     const std::vector<std::string>& specialSettings() override { return fSpecialSettings; }
     HepMC::GenEvent* apply(HepMC::GenEvent*) override;
+    HepMC3::GenEvent* apply(HepMC3::GenEvent*) override;
     void configureOnlyFor(int) override;
     void avoidTauLeptonicDecays() override {
       fAvoidTauLeptonicDecays = true;

GeneratorInterface/PhotosInterface/plugins/BuildFile.xml

@@ -4,6 +4,7 @@
   <use name="GeneratorInterface/PhotosInterface"/>
   <use name="photospp"/>
   <use name="hepmc"/>
+  <use name="hepmc3"/>
   <use name="clhep"/>
   <flags EDM_PLUGIN="1"/>
 </library>

GeneratorInterface/PhotosInterface/plugins/PhotosppInterface.cc

@@ -16,6 +16,7 @@ using namespace std;
 
 #include "Photos/Photos.h"
 #include "Photos/PhotosHepMCEvent.h"
+#include "Photos/PhotosHepMC3Event.h"
 
 CLHEP::HepRandomEngine* PhotosppInterface::fRandomEngine = nullptr;
 
@@ -285,6 +286,15 @@ HepMC::GenEvent* PhotosppInterface::apply(HepMC::GenEvent* evt) {
   return evt;
 }
 
+HepMC3::GenEvent* PhotosppInterface::apply(HepMC3::GenEvent* evt) {
+  Photospp::Photos::setRandomGenerator(PhotosppInterface::flat);
+  if (!fIsInitialized)
+    return evt;
+  Photospp::PhotosHepMC3Event PhotosEvt(evt);
+  PhotosEvt.process();
+  return evt;
+}
+
 double PhotosppInterface::flat() {
   if (!fRandomEngine) {
     throw cms::Exception("LogicError")

GeneratorInterface/Pythia8Interface/plugins/Pythia8HepMC3Hadronizer.cc

@@ -895,56 +895,72 @@ bool Pythia8HepMC3Hadronizer::residualDecay() {
   int NPartsBeforeDecays = pythiaEvent->size() - 1;  // do NOT count the very 1st "system" particle
                                                      // in Pythia8::Event record; it does NOT even
                                                      // get translated by the HepMCInterface to the
-                                                     // HepMC::GenEvent record !!!
-  //int NPartsAfterDecays = event().get()->particles_size();
-  int NPartsAfterDecays = 0;
-  for (auto p : (event3().get())->particles()) {
-    NPartsAfterDecays++;
-  }
+                                                     // HepMC3::GenEvent record !!!
 
+  int NPartsAfterDecays = ((event3().get())->particles()).size();
   if (NPartsAfterDecays == NPartsBeforeDecays)
     return true;
 
   bool result = true;
 
-// this below is to be rewritten in future development, then it will be uncommented
-#if 0
-  for (int ipart = NPartsAfterDecays; ipart > NPartsBeforeDecays; ipart--) {
-    HepMC::GenParticle *part = event().get()->barcode_to_particle(ipart);
-
-    if (part->status() == 1 && (fDecayer->particleData).canDecay(part->pdg_id())) {
-      fDecayer->event.reset();
-      Particle py8part(part->pdg_id(),
-                       93,
-                       0,
-                       0,
-                       0,
-                       0,
-                       0,
-                       0,
-                       part->momentum().x(),
-                       part->momentum().y(),
-                       part->momentum().z(),
-                       part->momentum().t(),
-                       part->generated_mass());
-      HepMC::GenVertex *ProdVtx = part->production_vertex();
-      py8part.vProd(ProdVtx->position().x(), ProdVtx->position().y(), ProdVtx->position().z(), ProdVtx->position().t());
-      py8part.tau((fDecayer->particleData).tau0(part->pdg_id()));
-      fDecayer->event.append(py8part);
-      int nentries = fDecayer->event.size();
-      if (!fDecayer->event[nentries - 1].mayDecay())
-        continue;
-      fDecayer->next();
-      int nentries1 = fDecayer->event.size();
-      if (nentries1 <= nentries)
-        continue;  //same number of particles, no decays...
-
-      part->set_status(2);
-
-      result = toHepMC.fill_next_event(*(fDecayer.get()), event().get(), -1, true, part);
+  for (const auto &p : (event3().get())->particles()) {
+    if (p->id() > NPartsBeforeDecays) {
+      if (p->status() == 1 && (fDecayer->particleData).canDecay(p->pid())) {
+        fDecayer->event.reset();
+        Particle py8part(p->pid(),
+                         93,
+                         0,
+                         0,
+                         0,
+                         0,
+                         0,
+                         0,
+                         p->momentum().x(),
+                         p->momentum().y(),
+                         p->momentum().z(),
+                         p->momentum().t(),
+                         p->generated_mass());
+
+        py8part.vProd(p->production_vertex()->position().x(),
+                      p->production_vertex()->position().y(),
+                      p->production_vertex()->position().z(),
+                      p->production_vertex()->position().t());
+
+        py8part.tau((fDecayer->particleData).tau0(p->pid()));
+        fDecayer->event.append(py8part);
+        int nentries = fDecayer->event.size();
+        if (!fDecayer->event[nentries - 1].mayDecay())
+          continue;
+        result = fDecayer->next();
+        int nentries1 = fDecayer->event.size();
+        if (nentries1 <= nentries)
+          continue;  //same number of particles, no decays...
+
+        p->set_status(2);
+
+        HepMC3::GenVertexPtr prod_vtx0 = make_shared<HepMC3::GenVertex>(  // neglect particle path to decay
+            HepMC3::FourVector(p->production_vertex()->position().x(),
+                               p->production_vertex()->position().y(),
+                               p->production_vertex()->position().z(),
+                               p->production_vertex()->position().t()));
+        prod_vtx0->add_particle_in(p);
+        (event3().get())->add_vertex(prod_vtx0);
+        HepMC3::GenParticle *pnew;
+        Pythia8::Event pyev = fDecayer->event;
+        double momFac = 1.;
+        for (int i = 2; i < pyev.size(); ++i) {
+          // Fill the particle.
+          pnew = new HepMC3::GenParticle(
+              HepMC3::FourVector(
+                  momFac * pyev[i].px(), momFac * pyev[i].py(), momFac * pyev[i].pz(), momFac * pyev[i].e()),
+              pyev[i].id(),
+              pyev[i].statusHepMC());
+          pnew->set_generated_mass(momFac * pyev[i].m());
+          prod_vtx0->add_particle_out(pnew);
+        }
+      }
     }
   }
-#endif
 
   return result;
 }