HepMC3 output for HadronizerFilter

mseidel42 · mseidel42 · commit ca88de935bd9 · 2025-03-12T13:10:12.000+01:00
diff --git a/GeneratorInterface/Core/interface/BaseHepMC3Filter.h b/GeneratorInterface/Core/interface/BaseHepMC3Filter.h
@@ -0,0 +1,22 @@
+#ifndef BaseHepMC3Filter_H
+#define BaseHepMC3Filter_H
+
+// base class for simple filter to run inside of HadronizerFilter for
+// multiple hadronization attempts on lhe events
+
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "SimDataFormats/GeneratorProducts/interface/HepMC3Product.h"
+
+//
+// class declaration
+//
+
+class BaseHepMC3Filter {
+public:
+  BaseHepMC3Filter();
+  ~BaseHepMC3Filter();
+
+  bool filter(const HepMC3::GenEvent* evt) { return true; };
+};
+
+#endif
diff --git a/GeneratorInterface/Core/interface/HadronizerFilter.h b/GeneratorInterface/Core/interface/HadronizerFilter.h
@@ -34,6 +34,7 @@
 // #include "GeneratorInterface/ExternalDecays/interface/ExternalDecayDriver.h"
 
 #include "GeneratorInterface/Core/interface/HepMCFilterDriver.h"
+#include "GeneratorInterface/Core/interface/HepMC3FilterDriver.h"
 
 // LHE Run
 #include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h"
@@ -44,10 +45,12 @@
 #include "GeneratorInterface/LHEInterface/interface/LHEEvent.h"
 
 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
+#include "SimDataFormats/GeneratorProducts/interface/HepMC3Product.h"
 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
 #include "SimDataFormats/GeneratorProducts/interface/GenLumiInfoHeader.h"
 #include "SimDataFormats/GeneratorProducts/interface/GenLumiInfoProduct.h"
 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
+#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct3.h"
 
 namespace edm {
   template <class HAD, class DEC>
@@ -81,11 +84,13 @@ namespace edm {
     // gen::ExternalDecayDriver* decayer_;
     Decayer* decayer_;
     HepMCFilterDriver* filter_;
+    HepMC3FilterDriver* filter3_;
     InputTag runInfoProductTag_;
     EDGetTokenT<LHERunInfoProduct> runInfoProductToken_;
     EDGetTokenT<LHEEventProduct> eventProductToken_;
     unsigned int counterRunInfoProducts_;
     unsigned int nAttempts_;
+    unsigned int ivhepmc = 2;
   };
 
   //------------------------------------------------------------------------
@@ -98,6 +103,7 @@ namespace edm {
         hadronizer_(ps),
         decayer_(nullptr),
         filter_(nullptr),
+        filter3_(nullptr),
         runInfoProductTag_(),
         runInfoProductToken_(),
         eventProductToken_(),
@@ -143,6 +149,7 @@ namespace edm {
     if (ps.exists("HepMCFilter")) {
       ParameterSet psfilter = ps.getParameter<ParameterSet>("HepMCFilter");
       filter_ = new HepMCFilterDriver(psfilter);
+      filter3_ = new HepMC3FilterDriver(psfilter);
     }
 
     //initialize setting for multiple hadronization attempts
@@ -156,8 +163,14 @@ namespace edm {
       usesResource(edm::uniqueSharedResourceName());
     }
 
-    produces<edm::HepMCProduct>("unsmeared");
-    produces<GenEventInfoProduct>();
+    ivhepmc = hadronizer_.getVHepMC();
+    if (ivhepmc == 2) {
+      produces<edm::HepMCProduct>("unsmeared");
+      produces<GenEventInfoProduct>();
+    } else if (ivhepmc == 3) {
+      produces<edm::HepMC3Product>("unsmeared");
+      produces<GenEventInfoProduct3>();
+    }
     produces<GenLumiInfoHeader, edm::Transition::BeginLuminosityBlock>();
     produces<GenLumiInfoProduct, edm::Transition::EndLuminosityBlock>();
     produces<GenRunInfoProduct, edm::Transition::EndRun>();
@@ -171,6 +184,8 @@ namespace edm {
       delete decayer_;
     if (filter_)
       delete filter_;
+    if (filter3_)
+      delete filter3_;
   }
 
   template <class HAD, class DEC>
@@ -186,7 +201,9 @@ namespace edm {
     ev.getByToken(eventProductToken_, product);
 
     std::unique_ptr<HepMC::GenEvent> finalEvent;
+    std::unique_ptr<HepMC3::GenEvent> finalEvent3;
     std::unique_ptr<GenEventInfoProduct> finalGenEventInfo;
+    std::unique_ptr<GenEventInfoProduct3> finalGenEventInfo3;
 
     //number of accepted events
     unsigned int naccept = 0;
@@ -208,8 +225,11 @@ namespace edm {
         continue;
 
       std::unique_ptr<HepMC::GenEvent> event(hadronizer_.getGenEvent());
-      if (!event.get())
-        continue;
+      std::unique_ptr<HepMC3::GenEvent> event3(hadronizer_.getGenEvent3());
+      if (ivhepmc == 2 && !event.get())
+        return false;
+      if (ivhepmc == 3 && !event3.get())
+        return false;
 
       // The external decay driver is being added to the system,
       // it should be called here
@@ -223,60 +243,105 @@ namespace edm {
         hadronizer_.setLHEEvent(std::move(lheEvent));
       }
 
-      if (!event.get())
-        continue;
+      if (ivhepmc == 2 && !event.get())
+        return false;
+      if (ivhepmc == 3 && !event3.get())
+        return false;
 
       // check and perform if there're any unstable particles after
       // running external decay packges
       //
       hadronizer_.resetEvent(std::move(event));
+      hadronizer_.resetEvent3(std::move(event3));
       if (!hadronizer_.residualDecay())
         continue;
 
       hadronizer_.finalizeEvent();
 
       event = hadronizer_.getGenEvent();
-      if (!event.get())
-        continue;
+      event3 = hadronizer_.getGenEvent3();
+      if (ivhepmc == 2 && !event.get())
+        return false;
+      if (ivhepmc == 3 && !event3.get())
+        return false;
+
+      if (ivhepmc == 2) {  // HepMC
+        event->set_event_number(ev.id().event());
+        std::unique_ptr<GenEventInfoProduct> genEventInfo(hadronizer_.getGenEventInfo());
+        if (!genEventInfo.get()) {
+          // create GenEventInfoProduct from HepMC event in case hadronizer didn't provide one
+          genEventInfo = std::make_unique<GenEventInfoProduct>(event.get());
+        }
 
-      event->set_event_number(ev.id().event());
+        //if HepMCFilter was specified, test event
+        if (filter_ && !filter_->filter(event.get(), genEventInfo->weight()))
+          continue;
 
-      std::unique_ptr<GenEventInfoProduct> genEventInfo(hadronizer_.getGenEventInfo());
-      if (!genEventInfo.get()) {
-        // create GenEventInfoProduct from HepMC event in case hadronizer didn't provide one
-        genEventInfo = std::make_unique<GenEventInfoProduct>(event.get());
-      }
+        ++naccept;
 
-      //if HepMCFilter was specified, test event
-      if (filter_ && !filter_->filter(event.get(), genEventInfo->weight()))
-        continue;
+        //keep the LAST accepted event (which is equivalent to choosing randomly from the accepted events)
+        finalEvent = std::move(event);
+        finalGenEventInfo = std::move(genEventInfo);
 
-      ++naccept;
+        if (!naccept)
+          return false;
 
-      //keep the LAST accepted event (which is equivalent to choosing randomly from the accepted events)
-      finalEvent = std::move(event);
-      finalGenEventInfo = std::move(genEventInfo);
-    }
+        //adjust event weights if necessary (in case input event was attempted multiple times)
+        if (nAttempts_ > 1) {
+          double multihadweight = double(naccept) / double(nAttempts_);
 
-    if (!naccept)
-      return false;
+          //adjust weight for GenEventInfoProduct
+          finalGenEventInfo->weights()[0] *= multihadweight;
 
-    //adjust event weights if necessary (in case input event was attempted multiple times)
-    if (nAttempts_ > 1) {
-      double multihadweight = double(naccept) / double(nAttempts_);
+          //adjust weight for HepMC GenEvent (used e.g for RIVET)
+          finalEvent->weights()[0] *= multihadweight;
+        }
 
-      //adjust weight for GenEventInfoProduct
-      finalGenEventInfo->weights()[0] *= multihadweight;
+        ev.put(std::move(finalGenEventInfo));
 
-      //adjust weight for HepMC GenEvent (used e.g for RIVET)
-      finalEvent->weights()[0] *= multihadweight;
-    }
+        std::unique_ptr<HepMCProduct> bare_product(new HepMCProduct());
+        bare_product->addHepMCData(finalEvent.release());
+        ev.put(std::move(bare_product), "unsmeared");
+
+      } else if (ivhepmc == 3) {  // HepMC3
+        event3->set_event_number(ev.id().event());
+        std::unique_ptr<GenEventInfoProduct3> genEventInfo3(hadronizer_.getGenEventInfo3());
+        if (!genEventInfo3.get()) {
+          // create GenEventInfoProduct3 from HepMC3 event in case hadronizer didn't provide one
+          genEventInfo3 = std::make_unique<GenEventInfoProduct3>(event3.get());
+        }
+
+        //if HepMCFilter was specified, test event
+        if (filter3_ && !filter3_->filter(event3.get(), genEventInfo3->weight()))
+          continue;
 
-    ev.put(std::move(finalGenEventInfo));
+        ++naccept;
 
-    std::unique_ptr<HepMCProduct> bare_product(new HepMCProduct());
-    bare_product->addHepMCData(finalEvent.release());
-    ev.put(std::move(bare_product), "unsmeared");
+        //keep the LAST accepted event (which is equivalent to choosing randomly from the accepted events)
+        finalEvent3 = std::move(event3);
+        finalGenEventInfo3 = std::move(genEventInfo3);
+
+        if (!naccept)
+          return false;
+
+        //adjust event weights if necessary (in case input event was attempted multiple times)
+        if (nAttempts_ > 1) {
+          double multihadweight = double(naccept) / double(nAttempts_);
+
+          //adjust weight for GenEventInfoProduct
+          finalGenEventInfo3->weights()[0] *= multihadweight;
+
+          //adjust weight for HepMC GenEvent (used e.g for RIVET)
+          finalEvent3->weights()[0] *= multihadweight;
+        }
+
+        ev.put(std::move(finalGenEventInfo3));
+
+        std::unique_ptr<HepMC3Product> bare_product(new HepMC3Product());
+        bare_product->addHepMCData(finalEvent3.release());
+        ev.put(std::move(bare_product), "unsmeared");
+      }
+    }
 
     return true;
   }
@@ -327,6 +392,8 @@ namespace edm {
       decayer_->statistics();
     if (filter_)
       filter_->statistics();
+    if (filter3_)
+      filter3_->statistics();
     lheRunInfo->statistics();
 
     std::unique_ptr<GenRunInfoProduct> griproduct(new GenRunInfoProduct(genRunInfo));
@@ -363,6 +430,9 @@ namespace edm {
     if (filter_) {
       filter_->resetStatistics();
     }
+    if (filter3_) {
+      filter3_->resetStatistics();
+    }
 
     if (!hadronizer_.initializeForExternalPartons())
       throw edm::Exception(errors::Configuration)
@@ -416,6 +486,17 @@ namespace edm {
                                                                    filter_->sumtotal_w2()));
       lumi.put(std::move(thisProduct));
     }
+    if (filter3_) {
+      std::unique_ptr<GenFilterInfo> thisProduct(new GenFilterInfo(filter3_->numEventsPassPos(),
+                                                                   filter3_->numEventsPassNeg(),
+                                                                   filter3_->numEventsTotalPos(),
+                                                                   filter3_->numEventsTotalNeg(),
+                                                                   filter3_->sumpass_w(),
+                                                                   filter3_->sumpass_w2(),
+                                                                   filter3_->sumtotal_w(),
+                                                                   filter3_->sumtotal_w2()));
+      lumi.put(std::move(thisProduct));
+    }
   }
 
 }  // namespace edm
diff --git a/GeneratorInterface/Core/interface/HepMC3FilterDriver.h b/GeneratorInterface/Core/interface/HepMC3FilterDriver.h
@@ -0,0 +1,45 @@
+#ifndef gen_HEPMC3FILTERDRIVER_H
+#define gen_HEPMC3FILTERDRIVER_H
+
+//class to select a HepMC3Filter to run with multiple hadronization attempts
+//inside HadronizerFilter
+
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "SimDataFormats/GeneratorProducts/interface/HepMC3Product.h"
+#include "GeneratorInterface/Core/interface/BaseHepMC3Filter.h"
+#include "SimDataFormats/GeneratorProducts/interface/GenFilterInfo.h"
+
+//
+// class declaration
+//
+
+class HepMC3FilterDriver {
+public:
+  HepMC3FilterDriver(const edm::ParameterSet&);
+  ~HepMC3FilterDriver();
+
+  bool filter(const HepMC3::GenEvent* evt, double weight = 1.);
+  void statistics() const;
+  void resetStatistics();
+
+  unsigned int numEventsPassPos() const { return numEventsPassPos_; }
+  unsigned int numEventsPassNeg() const { return numEventsPassNeg_; }
+  unsigned int numEventsTotalPos() const { return numEventsTotalPos_; }
+  unsigned int numEventsTotalNeg() const { return numEventsTotalNeg_; }
+  double sumpass_w() const { return sumpass_w_; }
+  double sumpass_w2() const { return sumpass_w2_; }
+  double sumtotal_w() const { return sumtotal_w_; }
+  double sumtotal_w2() const { return sumtotal_w2_; }
+
+private:
+  BaseHepMC3Filter* filter_;
+  unsigned int numEventsPassPos_;
+  unsigned int numEventsPassNeg_;
+  unsigned int numEventsTotalPos_;
+  unsigned int numEventsTotalNeg_;
+  double sumpass_w_;
+  double sumpass_w2_;
+  double sumtotal_w_;
+  double sumtotal_w2_;
+};
+#endif
diff --git a/GeneratorInterface/Core/src/BaseHepMC3Filter.cc b/GeneratorInterface/Core/src/BaseHepMC3Filter.cc
@@ -0,0 +1,5 @@
+#include "GeneratorInterface/Core/interface/BaseHepMC3Filter.h"
+
+BaseHepMC3Filter::BaseHepMC3Filter() {}
+
+BaseHepMC3Filter::~BaseHepMC3Filter() {}
diff --git a/GeneratorInterface/Core/src/HepMC3FilterDriver.cc b/GeneratorInterface/Core/src/HepMC3FilterDriver.cc
diff --git a/GeneratorInterface/Pythia8Interface/test/pythia8hmc3ex3_cfg.py b/GeneratorInterface/Pythia8Interface/test/pythia8hmc3ex3_cfg.py
