Add generator hepmc3 products

SimDataFormats/GeneratorProducts/BuildFile.xml

@@ -2,6 +2,7 @@
 <use name="FWCore/MessageLogger"/>
 <use name="DataFormats/Common"/>
 <use name="hepmc"/>
+<use name="hepmc3"/>
 <use name="xz"/>
 <use name="roothistmatrix"/>
 <export>

SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct3.h

@@ -0,0 +1,103 @@
+#ifndef SimDataFormats_GeneratorProducts_GenEventInfoProduct3_h
+#define SimDataFormats_GeneratorProducts_GenEventInfoProduct3_h
+
+#include <vector>
+#include <memory>
+
+#include "SimDataFormats/GeneratorProducts/interface/PdfInfo.h"
+
+namespace HepMC3 {
+  class GenEvent;
+}  // namespace HepMC3
+
+/** \class GenEventInfoProduct3
+ *
+ */
+
+class GenEventInfoProduct3 {
+public:
+  GenEventInfoProduct3();
+  GenEventInfoProduct3(const HepMC3::GenEvent *evt);
+  GenEventInfoProduct3(const GenEventInfoProduct3 &other);
+  GenEventInfoProduct3(GenEventInfoProduct3 &&other);
+  virtual ~GenEventInfoProduct3();
+
+  GenEventInfoProduct3 &operator=(const GenEventInfoProduct3 &other);
+  GenEventInfoProduct3 &operator=(GenEventInfoProduct3 &&other);
+
+  typedef gen::PdfInfo PDF;
+
+  // getters
+
+  std::vector<double> &weights() { return weights_; }
+  const std::vector<double> &weights() const { return weights_; }
+
+  double weight() const { return weights_.empty() ? 1.0 : weights_[0]; }
+
+  double weightProduct() const;
+
+  unsigned int signalProcessID() const { return signalProcessID_; }
+
+  double qScale() const { return qScale_; }
+  double alphaQCD() const { return alphaQCD_; }
+  double alphaQED() const { return alphaQED_; }
+
+  const PDF *pdf() const { return pdf_.get(); }
+  bool hasPDF() const { return pdf() != nullptr; }
+
+  const std::vector<double> &binningValues() const { return binningValues_; }
+  bool hasBinningValues() const { return !binningValues_.empty(); }
+
+  const std::vector<float> &DJRValues() const { return DJRValues_; }
+  bool hasDJRValues() const { return !DJRValues_.empty(); }
+
+  int nMEPartons() const { return nMEPartons_; }
+
+  int nMEPartonsFiltered() const { return nMEPartonsFiltered_; }
+
+  // setters
+
+  void setWeights(const std::vector<double> &weights) { weights_ = weights; }
+
+  void setSignalProcessID(unsigned int procID) { signalProcessID_ = procID; }
+
+  void setScales(double q = -1., double qcd = -1., double qed = -1.) { qScale_ = q, alphaQCD_ = qcd, alphaQED_ = qed; }
+
+  void setPDF(const PDF *pdf) { pdf_.reset(pdf ? new PDF(*pdf) : nullptr); }
+
+  void setBinningValues(const std::vector<double> &values) { binningValues_ = values; }
+
+  void setDJR(const std::vector<float> &values) { DJRValues_ = values; }
+
+  void setNMEPartons(int n) { nMEPartons_ = n; }
+
+  void setNMEPartonsFiltered(int n) { nMEPartonsFiltered_ = n; }
+
+private:
+  // HepMC3::GenEvent provides a list of weights
+  std::vector<double> weights_;
+
+  // generator-dependent process ID
+  unsigned int signalProcessID_;
+
+  // information about scales
+  double qScale_;
+  double alphaQCD_, alphaQED_;
+
+  // optional PDF info
+  std::unique_ptr<PDF> pdf_;
+
+  // If event was produced in bis, this contains
+  // the values that were used to define which
+  // bin the event belongs in
+  // This replaces the genEventScale, which only
+  // corresponds to Pythia pthat.  The RunInfo
+  // will contain the information what physical
+  // quantity these values actually belong to
+  std::vector<double> binningValues_;
+  std::vector<float> DJRValues_;
+  int nMEPartons_;
+  int nMEPartonsFiltered_;
+};
+
+#endif  // SimDataFormats_GeneratorProducts_GenEventInfoProduct3_h

SimDataFormats/GeneratorProducts/interface/HepMC3Product.h

@@ -0,0 +1,93 @@
+#ifndef SimDataFormats_GeneratorProducts_HepMC3Product_h
+#define SimDataFormats_GeneratorProducts_HepMC3Product_h
+
+/** \class HepMC3Product
+ *
+ *  \author Joanna Weng, Filip Moortgat, Mikhail Kirsanov
+ */
+
+#include "DataFormats/Common/interface/Ref.h"
+#include <TMatrixD.h>
+#include <HepMC3/GenEvent.h>
+#include <cstddef>
+
+namespace HepMC3 {
+  class FourVector;
+  class GenParticle;
+  class GenVertex;
+}  // namespace HepMC3
+
+namespace edm {
+  class HepMC3Product {
+  public:
+    HepMC3Product() : evt_(nullptr), isVtxGenApplied_(false), isVtxBoostApplied_(false), isPBoostApplied_(false) {}
+
+    explicit HepMC3Product(HepMC3::GenEvent *evt);
+    virtual ~HepMC3Product();
+
+    void addHepMCData(HepMC3::GenEvent *evt);
+
+    void applyVtxGen(HepMC3::FourVector const *vtxShift) { applyVtxGen(*vtxShift); }
+    void applyVtxGen(HepMC3::FourVector const &vtxShift);
+
+    void boostToLab(TMatrixD const *lorentz, std::string const &type);
+
+    const HepMC3::GenEvent &getHepMCData() const;
+
+    const HepMC3::GenEvent *GetEvent() const { return evt_; }
+
+    bool isVtxGenApplied() const { return isVtxGenApplied_; }
+    bool isVtxBoostApplied() const { return isVtxBoostApplied_; }
+    bool isPBoostApplied() const { return isPBoostApplied_; }
+
+    HepMC3Product(HepMC3Product const &orig);
+    HepMC3Product &operator=(HepMC3Product const &other);
+    HepMC3Product(HepMC3Product &&orig);
+    HepMC3Product &operator=(HepMC3Product &&other);
+    void swap(HepMC3Product &other);
+
+  private:
+    HepMC3::GenEvent *evt_;
+
+    bool isVtxGenApplied_;
+    bool isVtxBoostApplied_;
+    bool isPBoostApplied_;
+  };
+
+  // This allows edm::Refs to work with HepMC3Product
+  namespace refhelper {
+    template <>
+    struct FindTrait<edm::HepMC3Product, HepMC3::GenParticle> {
+      struct Find {
+        using first_argument_type = edm::HepMC3Product const &;
+        using second_argument_type = int;
+        using result_type = HepMC3::GenParticle const *;
+
+        result_type operator()(first_argument_type iContainer, second_argument_type iBarCode) {
+          //return iContainer.getHepMCData().barcode_to_particle(iBarCode);
+          return nullptr;
+        }
+      };
+
+      typedef Find value;
+    };
+
+    template <>
+    struct FindTrait<edm::HepMC3Product, HepMC3::GenVertex> {
+      struct Find {
+        using first_argument_type = edm::HepMC3Product const &;
+        using second_argument_type = int;
+        using result_type = HepMC3::GenVertex const *;
+
+        result_type operator()(first_argument_type iContainer, second_argument_type iBarCode) {
+          //return iContainer.getHepMCData().barcode_to_vertex(iBarCode);
+          return nullptr;
+        }
+      };
+
+      typedef Find value;
+    };
+  }  // namespace refhelper
+}  // namespace edm
+
+#endif  // SimDataFormats_GeneratorProducts_HepMC3Product_h

SimDataFormats/GeneratorProducts/src/GenEventInfoProduct3.cc

@@ -0,0 +1,103 @@
+#include <functional>
+#include <numeric>
+using std::ptrdiff_t;
+
+#include <HepMC3/GenEvent.h>
+// #include <HepMC3/WeightContainer.h>
+#include <HepMC3/GenPdfInfo.h>
+
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+
+#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct3.h"
+
+using namespace edm;
+using namespace std;
+
+GenEventInfoProduct3::GenEventInfoProduct3()
+    : signalProcessID_(0), qScale_(-1.), alphaQCD_(-1.), alphaQED_(-1.), nMEPartons_(-1), nMEPartonsFiltered_(-1) {}
+
+GenEventInfoProduct3::GenEventInfoProduct3(const HepMC3::GenEvent *evt)
+    : weights_(evt->weights().begin(), evt->weights().end()), nMEPartons_(-1), nMEPartonsFiltered_(-1) {
+  std::shared_ptr<HepMC3::IntAttribute> A_signal_process_id = evt->attribute<HepMC3::IntAttribute>("signal_process_id");
+  std::shared_ptr<HepMC3::DoubleAttribute> A_event_scale = evt->attribute<HepMC3::DoubleAttribute>("event_scale");
+  std::shared_ptr<HepMC3::DoubleAttribute> A_alphaQCD = evt->attribute<HepMC3::DoubleAttribute>("alphaQCD");
+  std::shared_ptr<HepMC3::DoubleAttribute> A_alphaQED = evt->attribute<HepMC3::DoubleAttribute>("alphaQED");
+  //std::shared_ptr<HepMC3::IntAttribute> A_mpi = evt->attribute<HepMC3::IntAttribute>("mpi");
+
+  signalProcessID_ = A_signal_process_id ? (A_signal_process_id->value()) : 0;
+  qScale_ = A_event_scale ? (A_event_scale->value()) : 0.0;
+  alphaQCD_ = A_alphaQCD ? (A_alphaQCD->value()) : 0.0;
+  alphaQED_ = A_alphaQED ? (A_alphaQED->value()) : 0.0;
+
+  std::shared_ptr<HepMC3::GenPdfInfo> A_pdf = evt->attribute<HepMC3::GenPdfInfo>("GenPdfInfo");
+  if (A_pdf) {
+    PDF pdf;
+    pdf.id = std::make_pair(A_pdf->parton_id[0], A_pdf->parton_id[1]);
+    pdf.x = std::make_pair(A_pdf->x[0], A_pdf->x[1]);
+    pdf.xPDF = std::make_pair(A_pdf->xf[0], A_pdf->xf[1]);
+    pdf.scalePDF = A_pdf->scale;
+    setPDF(&pdf);
+  }
+}
+
+GenEventInfoProduct3::GenEventInfoProduct3(GenEventInfoProduct3 const &other)
+    : weights_(other.weights_),
+      signalProcessID_(other.signalProcessID_),
+      qScale_(other.qScale_),
+      alphaQCD_(other.alphaQCD_),
+      alphaQED_(other.alphaQED_),
+      binningValues_(other.binningValues_),
+      DJRValues_(other.DJRValues_),
+      nMEPartons_(other.nMEPartons_),
+      nMEPartonsFiltered_(other.nMEPartons_) {
+  setPDF(other.pdf());
+}
+
+GenEventInfoProduct3::GenEventInfoProduct3(GenEventInfoProduct3 &&other)
+    : weights_(std::move(other.weights_)),
+      signalProcessID_(other.signalProcessID_),
+      qScale_(other.qScale_),
+      alphaQCD_(other.alphaQCD_),
+      alphaQED_(other.alphaQED_),
+      pdf_(other.pdf_.release()),
+      binningValues_(std::move(other.binningValues_)),
+      DJRValues_(std::move(other.DJRValues_)),
+      nMEPartons_(other.nMEPartons_),
+      nMEPartonsFiltered_(other.nMEPartons_) {}
+
+GenEventInfoProduct3::~GenEventInfoProduct3() {}
+
+GenEventInfoProduct3 &GenEventInfoProduct3::operator=(GenEventInfoProduct3 const &other) {
+  weights_ = other.weights_;
+  signalProcessID_ = other.signalProcessID_;
+  qScale_ = other.qScale_;
+  alphaQCD_ = other.alphaQCD_;
+  alphaQED_ = other.alphaQED_;
+  binningValues_ = other.binningValues_;
+  DJRValues_ = other.DJRValues_;
+  nMEPartons_ = other.nMEPartons_;
+  nMEPartonsFiltered_ = other.nMEPartonsFiltered_;
+
+  setPDF(other.pdf());
+
+  return *this;
+}
+
+GenEventInfoProduct3 &GenEventInfoProduct3::operator=(GenEventInfoProduct3 &&other) {
+  weights_ = std::move(other.weights_);
+  signalProcessID_ = other.signalProcessID_;
+  qScale_ = other.qScale_;
+  alphaQCD_ = other.alphaQCD_;
+  alphaQED_ = other.alphaQED_;
+  binningValues_ = std::move(other.binningValues_);
+  DJRValues_ = std::move(other.DJRValues_);
+  nMEPartons_ = other.nMEPartons_;
+  nMEPartonsFiltered_ = other.nMEPartonsFiltered_;
+  pdf_ = std::move(other.pdf_);
+
+  return *this;
+}
+
+double GenEventInfoProduct3::weightProduct() const {
+  return std::accumulate(weights_.begin(), weights_.end(), 1., std::multiplies<double>());
+}