Hepmc3tog4, addind the codes that provide a possibility to transfer HepMC3Product from the GEN step to the SIM step

GeneratorInterface/Core/plugins/GeneratorSmearedProducer.cc

@@ -5,6 +5,7 @@
 #include "FWCore/ParameterSet/interface/ParameterSet.h"
 
 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
+#include "SimDataFormats/GeneratorProducts/interface/HepMC3Product.h"
 #include "FWCore/Framework/interface/global/EDProducer.h"
 #include "FWCore/Utilities/interface/EDGetToken.h"
 #include "FWCore/Framework/interface/MakerMacros.h"
@@ -16,7 +17,6 @@ namespace edm {
   class ConfigurationDescriptions;
   class Event;
   class EventSetup;
-  class HepMCProduct;
 }  // namespace edm
 
 class GeneratorSmearedProducer : public edm::global::EDProducer<> {
@@ -29,14 +29,19 @@ class GeneratorSmearedProducer : public edm::global::EDProducer<> {
 private:
   const edm::EDGetTokenT<edm::HepMCProduct> newToken_;
   const edm::EDGetTokenT<edm::HepMCProduct> oldToken_;
+  const edm::EDGetTokenT<edm::HepMC3Product> Token3_;
 };
 
 GeneratorSmearedProducer::GeneratorSmearedProducer(edm::ParameterSet const& ps)
     : newToken_(consumes<edm::HepMCProduct>(ps.getUntrackedParameter<edm::InputTag>("currentTag"))),
-      oldToken_(consumes<edm::HepMCProduct>(ps.getUntrackedParameter<edm::InputTag>("previousTag"))) {
+      oldToken_(consumes<edm::HepMCProduct>(ps.getUntrackedParameter<edm::InputTag>("previousTag"))),
+      Token3_(consumes<edm::HepMC3Product>(ps.getUntrackedParameter<edm::InputTag>("currentTag"))) {
   // This producer produces a HepMCProduct, a copy of the original one
-  // It is used for backwaerd compatibility
+  // It is used for backward compatibility
+  // If HepMC3Product exists, it produces its copy
+  // It adds "generatorSmeared" to description, which is needed for further processing
   produces<edm::HepMCProduct>();
+  produces<edm::HepMC3Product>();
 }
 
 void GeneratorSmearedProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& es) const {
@@ -49,6 +54,13 @@ void GeneratorSmearedProducer::produce(edm::StreamID, edm::Event& iEvent, const
     std::unique_ptr<edm::HepMCProduct> theCopy(new edm::HepMCProduct(*theHepMCProduct));
     iEvent.put(std::move(theCopy));
   }
+
+  edm::Handle<edm::HepMC3Product> theHepMC3Product;
+  found = iEvent.getByToken(Token3_, theHepMC3Product);
+  if (found) {
+    std::unique_ptr<edm::HepMC3Product> theCopy3(new edm::HepMC3Product(*theHepMC3Product));
+    iEvent.put(std::move(theCopy3));
+  }
 }
 
 void GeneratorSmearedProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {

GeneratorInterface/Pythia8Interface/test/pythia8test/runtest

@@ -4,4 +4,4 @@ rm -f p8test.html
 cmsRun p8test1step1_cfg.py
 cmsRun ZJetsTest_cfg.py
 cp testi.dat test.dat
-cmpr.sh test.dat html > p8test.html
+./cmpr.sh test.dat html > p8test.html

IOMC/EventVertexGenerators/interface/BaseEvtVtxGenerator.h

@@ -18,7 +18,8 @@ namespace CLHEP {
 
 namespace edm {
   class HepMCProduct;
-}
+  class HepMC3Product;
+}  // namespace edm
 
 class BaseEvtVtxGenerator : public edm::stream::EDProducer<> {
 public:
@@ -39,6 +40,7 @@ class BaseEvtVtxGenerator : public edm::stream::EDProducer<> {
 
 private:
   edm::EDGetTokenT<edm::HepMCProduct> sourceToken;
+  edm::EDGetTokenT<edm::HepMC3Product> sourceToken3;
 };
 
 #endif

IOMC/EventVertexGenerators/src/BaseEvtVtxGenerator.cc

@@ -8,7 +8,7 @@
 #include "FWCore/ParameterSet/interface/ParameterSet.h"
 
 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
-#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
+#include "SimDataFormats/GeneratorProducts/interface/HepMC3Product.h"
 
 #include "FWCore/ServiceRegistry/interface/Service.h"
 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
@@ -35,7 +35,9 @@ BaseEvtVtxGenerator::BaseEvtVtxGenerator(const ParameterSet& pset) {
                                              "in the configuration file or remove the modules that require it.";
   }
 
+  sourceToken3 = consumes<edm::HepMC3Product>(pset.getParameter<edm::InputTag>("src"));
   sourceToken = consumes<edm::HepMCProduct>(pset.getParameter<edm::InputTag>("src"));
+  produces<edm::HepMC3Product>();
   produces<edm::HepMCProduct>();
 }
 
@@ -47,20 +49,37 @@ void BaseEvtVtxGenerator::produce(Event& evt, const EventSetup&) {
 
   Handle<HepMCProduct> HepUnsmearedMCEvt;
 
-  evt.getByToken(sourceToken, HepUnsmearedMCEvt);
+  bool found = evt.getByToken(sourceToken, HepUnsmearedMCEvt);
+
+  if (found) {  // HepMC event exists
+
+    // Make a copy
+    HepMC::GenEvent* genevt = new HepMC::GenEvent(*HepUnsmearedMCEvt->GetEvent());
+
+    std::unique_ptr<edm::HepMCProduct> HepMCEvt(new edm::HepMCProduct(genevt));
+    // generate new vertex & apply the shift
+    //
+    HepMCEvt->applyVtxGen(newVertex(engine));
 
-  // Copy the HepMC::GenEvent
-  HepMC::GenEvent* genevt = new HepMC::GenEvent(*HepUnsmearedMCEvt->GetEvent());
-  std::unique_ptr<edm::HepMCProduct> HepMCEvt(new edm::HepMCProduct(genevt));
-  // generate new vertex & apply the shift
-  //
-  HepMCEvt->applyVtxGen(newVertex(engine));
+    HepMCEvt->boostToLab(GetInvLorentzBoost(), "vertex");
+    HepMCEvt->boostToLab(GetInvLorentzBoost(), "momentum");
 
-  //HepMCEvt->LorentzBoost( 0., 142.e-6 );
-  HepMCEvt->boostToLab(GetInvLorentzBoost(), "vertex");
-  HepMCEvt->boostToLab(GetInvLorentzBoost(), "momentum");
+    evt.put(std::move(HepMCEvt));
 
-  evt.put(std::move(HepMCEvt));
+  } else {  // no HepMC event, try to get HepMC3 event
+
+    Handle<HepMC3Product> HepUnsmearedMCEvt3;
+    found = evt.getByToken(sourceToken3, HepUnsmearedMCEvt3);
+
+    if (!found)
+      throw cms::Exception("ProductAbsent") << "No HepMCProduct, tried to get HepMC3Product, but it is also absent.";
+
+    HepMC3::GenEvent* genevt3 = new HepMC3::GenEvent();
+    genevt3->read_data(*HepUnsmearedMCEvt3->GetEvent());
+    HepMC3Product* productcopy3 = new HepMC3Product(genevt3);  // For the moment do not really smear HepMC3
+    std::unique_ptr<edm::HepMC3Product> HepMC3Evt(productcopy3);
+    evt.put(std::move(HepMC3Evt));
+  }
 
   return;
 }

PhysicsTools/HepMCCandAlgos/plugins/GenParticleProducer.cc

@@ -196,11 +196,15 @@ void GenParticleProducer::produce(StreamID, Event& evt, const EventSetup& es) co
     LogDebug("GenParticleProducer") << "totalSize : " << totalSize << endl;
   } else {
     Handle<HepMCProduct> mcp;
-    evt.getByToken(srcToken_, mcp);
-    mc = mcp->GetEvent();
-    if (mc == nullptr)
-      throw edm::Exception(edm::errors::InvalidReference) << "HepMC has null pointer to GenEvent" << endl;
-    totalSize = mc->particles_size();
+    bool found = evt.getByToken(srcToken_, mcp);
+    if (found) {
+      mc = mcp->GetEvent();
+      if (mc == nullptr)
+        throw edm::Exception(edm::errors::InvalidReference) << "HepMC has null pointer to GenEvent" << endl;
+      totalSize = mc->particles_size();
+    } else {
+      totalSize = 0;
+    }
   }
 
   // initialise containers
@@ -281,28 +285,30 @@ void GenParticleProducer::produce(StreamID, Event& evt, const EventSetup& es) co
       offset += num_particles;
     }
   } else {
-    auto origin = (*mc->vertices_begin())->position();
-    xyz0Ptr->SetXYZ(origin.x() * mmToCm, origin.y() * mmToCm, origin.z() * mmToCm);
-    *t0Ptr = origin.t() * mmToNs;
-    fillIndices(mc, particles, *barCodeVector, 0, barcodes);
-
-    // fill output collection and save association
-    for (size_t i = 0; i < particles.size(); ++i) {
-      const HepMC::GenParticle* part = particles[i];
-      reco::GenParticle& cand = cands[i];
-      // convert HepMC::GenParticle to new reco::GenParticle
-      convertParticle(cand, part, id2Charge);
-      cand.resetDaughters(ref.id());
-    }
+    if (totalSize) {
+      auto origin = (*mc->vertices_begin())->position();
+      xyz0Ptr->SetXYZ(origin.x() * mmToCm, origin.y() * mmToCm, origin.z() * mmToCm);
+      *t0Ptr = origin.t() * mmToNs;
+      fillIndices(mc, particles, *barCodeVector, 0, barcodes);
 
-    // fill references to daughters
-    for (size_t d = 0; d < cands.size(); ++d) {
-      const HepMC::GenParticle* part = particles[d];
-      const GenVertex* productionVertex = part->production_vertex();
-      // search barcode map and attach daughters
-      if (productionVertex != nullptr)
-        fillDaughters(cands, part, ref, d, barcodes);
-      cands[d].setCollisionId(0);
+      // fill output collection and save association
+      for (size_t i = 0; i < particles.size(); ++i) {
+        const HepMC::GenParticle* part = particles[i];
+        reco::GenParticle& cand = cands[i];
+        // convert HepMC::GenParticle to new reco::GenParticle
+        convertParticle(cand, part, id2Charge);
+        cand.resetDaughters(ref.id());
+      }
+
+      // fill references to daughters
+      for (size_t d = 0; d < cands.size(); ++d) {
+        const HepMC::GenParticle* part = particles[d];
+        const GenVertex* productionVertex = part->production_vertex();
+        // search barcode map and attach daughters
+        if (productionVertex != nullptr)
+          fillDaughters(cands, part, ref, d, barcodes);
+        cands[d].setCollisionId(0);
+      }
     }
   }
 
@@ -353,7 +359,6 @@ bool GenParticleProducer::fillDaughters(reco::GenParticleCollection& cands,
       cands[index].addMother(GenParticleRef(ref, m));
     }
   }
-
   return true;
 }
 

SimG4Core/Application/interface/RunManagerMTWorker.h

@@ -7,6 +7,7 @@
 #include "SimDataFormats/Forward/interface/LHCTransportLinkContainer.h"
 
 #include "SimG4Core/Generators/interface/Generator.h"
+#include "SimG4Core/Generators/interface/Generator3.h"
 #include "SimG4Core/Notification/interface/TmpSimEvent.h"
 
 #include "MagneticField/Engine/interface/MagneticField.h"
@@ -22,6 +23,7 @@ namespace edm {
   class EventSetup;
   class ConsumesCollector;
   class HepMCProduct;
+  class HepMC3Product;
 }  // namespace edm
 
 class Generator;
@@ -89,7 +91,9 @@ class RunManagerMTWorker {
   void DumpMagneticField(const G4Field*, const std::string&) const;
 
   Generator m_generator;
+  Generator3 m_generator3;
   edm::EDGetTokenT<edm::HepMCProduct> m_InToken;
+  edm::EDGetTokenT<edm::HepMC3Product> m_InToken3;
   edm::EDGetTokenT<edm::HepMCProduct> m_LHCToken;
   edm::EDGetTokenT<edm::LHCTransportLinkContainer> m_theLHCTlinkToken;
   edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> m_MagField;

SimG4Core/Application/src/RunManagerMTWorker.cc

@@ -35,7 +35,9 @@
 #include "SimG4Core/MagneticField/interface/CMSFieldManager.h"
 
 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
+#include "SimDataFormats/GeneratorProducts/interface/HepMC3Product.h"
 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
+#include "HepMC3/Print.h"
 
 #include "SimG4Core/Physics/interface/PhysicsList.h"
 
@@ -152,8 +154,11 @@ struct RunManagerMTWorker::TLSData {
 
 RunManagerMTWorker::RunManagerMTWorker(const edm::ParameterSet& p, edm::ConsumesCollector&& iC)
     : m_generator(p.getParameter<edm::ParameterSet>("Generator")),
+      m_generator3(p.getParameter<edm::ParameterSet>("Generator")),
       m_InToken(iC.consumes<edm::HepMCProduct>(
           p.getParameter<edm::ParameterSet>("Generator").getParameter<edm::InputTag>("HepMCProductLabel"))),
+      m_InToken3(iC.consumes<edm::HepMC3Product>(
+          p.getParameter<edm::ParameterSet>("Generator").getParameter<edm::InputTag>("HepMCProductLabel"))),
       m_theLHCTlinkToken(iC.consumes<edm::LHCTransportLinkContainer>(p.getParameter<edm::InputTag>("theLHCTlinkTag"))),
       m_nonBeam(p.getParameter<bool>("NonBeamEvent")),
       m_UseG4EventManager(p.getParameter<bool>("UseG4EventManager")),
@@ -572,6 +577,7 @@ TmpSimEvent* RunManagerMTWorker::produce(const edm::Event& inpevt,
   }
 
   // event and primary
+
   m_tls->currentEvent.reset(generateEvent(inpevt));
   m_simEvent.clear();
   m_simEvent.setHepEvent(m_generator.genEvent());
@@ -645,26 +651,48 @@ G4Event* RunManagerMTWorker::generateEvent(const edm::Event& inpevt) {
   G4int evtid = (G4int)inpevt.id().event();
   G4Event* evt = new G4Event(evtid);
 
-  edm::Handle<edm::HepMCProduct> HepMCEvt;
-  inpevt.getByToken(m_InToken, HepMCEvt);
-
-  m_generator.setGenEvent(HepMCEvt->GetEvent());
-
   // required to reset the GenParticle Id for particles transported
   // along the beam pipe to their original value for SimTrack creation
   resetGenParticleId(inpevt);
 
-  if (!m_nonBeam) {
-    m_generator.HepMC2G4(HepMCEvt->GetEvent(), evt);
-    if (m_LHCTransport) {
-      edm::Handle<edm::HepMCProduct> LHCMCEvt;
-      inpevt.getByToken(m_LHCToken, LHCMCEvt);
-      m_generator.nonCentralEvent2G4(LHCMCEvt->GetEvent(), evt);
+  edm::Handle<edm::HepMCProduct> HepMCEvt;
+  bool found = inpevt.getByToken(m_InToken, HepMCEvt);
+
+  if (found) {  // HepMC event exists
+
+    m_generator.setGenEvent(HepMCEvt->GetEvent());
+
+    if (!m_nonBeam) {
+      m_generator.HepMC2G4(HepMCEvt->GetEvent(), evt);
+      if (m_LHCTransport) {
+        edm::Handle<edm::HepMCProduct> LHCMCEvt;
+        inpevt.getByToken(m_LHCToken, LHCMCEvt);
+        m_generator.nonCentralEvent2G4(LHCMCEvt->GetEvent(), evt);
+      }
+    } else {
+      m_generator.nonCentralEvent2G4(HepMCEvt->GetEvent(), evt);
+    }
+
+  } else {  // no HepMC event, try to get HepMC3 event
+
+    edm::Handle<edm::HepMC3Product> HepMCEvt3;
+    inpevt.getByToken(m_InToken3, HepMCEvt3);
+
+    HepMC3::GenEvent* genevt3 = new HepMC3::GenEvent();
+    genevt3->read_data(*HepMCEvt3->GetEvent());
+    m_generator3.setGenEvent(genevt3);
+
+    if (!m_nonBeam) {
+      m_generator3.HepMC2G4(genevt3, evt);
+      if (m_LHCTransport) {
+        edm::Handle<edm::HepMC3Product> LHCMCEvt;
+        inpevt.getByToken(m_LHCToken, LHCMCEvt);
+        //m_generator3.nonCentralEvent2G4(LHCMCEvt->GetEvent(), evt);
+      }
+    } else {
+      //m_generator3.nonCentralEvent2G4(HepMCEvt->GetEvent(), evt);
     }
-  } else {
-    m_generator.nonCentralEvent2G4(HepMCEvt->GetEvent(), evt);
   }
-
   return evt;
 }
 

SimG4Core/Generators/BuildFile.xml

@@ -5,6 +5,7 @@
 <use name="FWCore/MessageLogger"/>
 <use name="boost"/>
 <use name="hepmc"/>
+<use name="hepmc3"/>
 <use name="heppdt"/>
 <use name="geant4core"/>
 <use name="rootmath"/>