HYDJET, PYQUEN and the new embedding

Configuration/Applications/python/ConfigBuilder.py

@@ -1379,6 +1379,8 @@ def prepare_GEN(self, sequence = None):
             if self._options.scenario == 'HeavyIons':
                 if self._options.pileup=='HiMixGEN':
                     self.loadAndRemember("Configuration/StandardSequences/GeneratorMix_cff")
+                elif self._options.pileup=='HiMixEmbGEN':
+                    self.loadAndRemember("Configuration/StandardSequences/GeneratorEmbMix_cff")
                 else:
                     self.loadAndRemember("Configuration/StandardSequences/GeneratorHI_cff")
 
@@ -2182,6 +2184,9 @@ def prepare(self, doChecking = False):
                 self.pythonCfgCode +='\n'
                 self.pythonCfgCode +=dumpPython(self.process,object)
 
+        if self._options.pileup=='HiMixEmbGEN':
+            self.pythonCfgCode += "\nprocess.generator.embeddingMode=cms.bool(True)\n"
+
         # dump all paths
         self.pythonCfgCode += "\n# Path and EndPath definitions\n"
         for path in self.process.paths:
@@ -2246,7 +2251,7 @@ def prepare(self, doChecking = False):
             MassReplaceInputTag(self.process, new="rawDataMapperByLabel", old="rawDataCollector")
 
         # special treatment in case of production filter sequence 2/2
-        if self.productionFilterSequence:
+        if self.productionFilterSequence and not (self._options.pileup=='HiMixEmbGEN'):
             self.pythonCfgCode +='# filter all path with the production filter sequence\n'
             self.pythonCfgCode +='for path in process.paths:\n'
             if len(self.conditionalPaths):

Configuration/Generator/python/Pyquen2013Settings_cff.py

@@ -28,7 +28,6 @@
 pyquenParameters  = cms.PSet(doIsospin = cms.bool(True),
                              angularSpectrumSelector = cms.int32(0), ## angular emitted gluon spectrum :
                              embeddingMode = cms.bool(False),
-                             backgroundLabel = cms.InputTag("generator","unsmeared") ## ineffective in no mixing
                              )
 
 hydjetParameters = cms.PSet(sigmaInelNN = cms.double(64),

Configuration/Generator/python/Pyquen2015Settings_cff.py

@@ -28,7 +28,6 @@
 pyquenParameters  = cms.PSet(doIsospin = cms.bool(True),
                              angularSpectrumSelector = cms.int32(0), ## angular emitted gluon spectrum
                              embeddingMode = cms.bool(False),
-                             backgroundLabel = cms.InputTag("generator","unsmeared") ## ineffective in no mixing
                              )
 
 hydjetParameters = cms.PSet(sigmaInelNN = cms.double(70),

Configuration/Generator/python/PyquenDefaultSettings_cff.py

@@ -26,7 +26,6 @@
 pyquenParameters  = cms.PSet(doIsospin = cms.bool(True),
                              angularSpectrumSelector = cms.int32(0), ## angular emitted gluon spectrum :
                              embeddingMode = cms.bool(False),
-                             backgroundLabel = cms.InputTag("generator","unsmeared") ## ineffective in no mixing
                              )
 
 hydjetParameters = cms.PSet(sigmaInelNN = cms.double(58),

Configuration/StandardSequences/python/GeneratorEmbMix_cff.py

@@ -0,0 +1,4 @@
+import FWCore.ParameterSet.Config as cms
+from Configuration.StandardSequences.GeneratorMix_cff import *
+
+pgen = cms.Sequence(cms.SequencePlaceholder("mix")+cms.SequencePlaceholder("generator")+cms.SequencePlaceholder("randomEngineStateProducer")+VertexSmearing+GenSmeared+genParticles+hiGenJets)

Configuration/StandardSequences/python/Mixing.py

@@ -53,6 +53,7 @@ def addMixingScenario(label,dict):
 addMixingScenario("E8TeV_AVE_10_BX_25ns_300ns_spread",{'file':'SimGeneral.MixingModule.mix_E8TeV_AVE_10_BX_25ns_300ns_spread_cfi'})
 addMixingScenario("HiMix",{'file': 'SimGeneral.MixingModule.HiMix_cff'})
 addMixingScenario("HiMixGEN",{'file': 'SimGeneral.MixingModule.HiMixGEN_cff'})
+addMixingScenario("HiMixEmbGEN",{'file': 'SimGeneral.MixingModule.HiMixEmbGEN_cff'})
 addMixingScenario("HiMixNoPU",{'file': 'SimGeneral.MixingModule.HiMixNoPU_cff'})
 addMixingScenario("HighLumiPileUp",{'file': 'SimGeneral.MixingModule.mixHighLumPU_cfi'})
 addMixingScenario("InitialLumiPileUp",{'file': 'SimGeneral.MixingModule.mixInitialLumPU_cfi'})

GeneratorInterface/Core/interface/GeneratorFilter.h

@@ -61,10 +61,13 @@ namespace edm {
     void preallocThreads(unsigned int iThreads) override;
 
   private:
+    // two-phase construction to allow specialization of the constructor
+    void init(ParameterSet const& ps);
+
     Hadronizer hadronizer_;
     //gen::ExternalDecayDriver* decayer_;
-    Decayer* decayer_;
-    unsigned int nEventsInLumiBlock_;
+    Decayer* decayer_ = nullptr;
+    unsigned int nEventsInLumiBlock_ = 0;
     unsigned int nThreads_{1};
     bool initialized_ = false;
   };
@@ -74,8 +77,12 @@ namespace edm {
   // Implementation
 
   template <class HAD, class DEC>
-  GeneratorFilter<HAD, DEC>::GeneratorFilter(ParameterSet const& ps)
-      : EDFilter(), hadronizer_(ps), decayer_(nullptr), nEventsInLumiBlock_(0) {
+  GeneratorFilter<HAD, DEC>::GeneratorFilter(ParameterSet const& ps) : hadronizer_(ps) {
+    init(ps);
+  }
+
+  template <class HAD, class DEC>
+  void GeneratorFilter<HAD, DEC>::init(ParameterSet const& ps) {
     // TODO:
     // Put the list of types produced by the filters here.
     // The current design calls for:

GeneratorInterface/HydjetInterface/BuildFile.xml

@@ -1,12 +1,12 @@
-<use name="boost"/>
-<use name="GeneratorInterface/Core"/>
-<use name="FWCore/Concurrency"/>
-<use name="FWCore/Framework"/>
-<use name="SimDataFormats/GeneratorProducts"/>
-<use name="GeneratorInterface/Pythia6Interface"/>
-<use name="GeneratorInterface/PyquenInterface"/>
-<use name="GeneratorInterface/ExternalDecays"/>
-<use name="f77compiler"/>
+<use   name="boost"/>
+<use   name="GeneratorInterface/Core"/>
+<use   name="FWCore/Concurrency"/>
+<use   name="FWCore/Framework"/>
+<use   name="SimDataFormats/GeneratorProducts"/>
+<use   name="GeneratorInterface/Pythia6Interface"/>
+<use   name="GeneratorInterface/ExternalDecays"/>
+<use   name="f77compiler"/>
+<use   name="hydjet"/>
 <export>
-  <lib name="1"/>
+  <lib   name="1"/>
 </export>

GeneratorInterface/HydjetInterface/interface/HydjetGeneratorFilter.h

@@ -5,6 +5,14 @@
 #include "GeneratorInterface/Core/interface/GeneratorFilter.h"
 #include "GeneratorInterface/ExternalDecays/interface/ExternalDecayDriver.h"
 
+namespace edm {
+  template <>
+  GeneratorFilter<gen::HydjetHadronizer, gen::ExternalDecayDriver>::GeneratorFilter(ParameterSet const& ps)
+      : hadronizer_(ps, consumesCollector()) {
+    init(ps);
+  }
+}  // namespace edm
+
 namespace gen {
   typedef edm::GeneratorFilter<gen::HydjetHadronizer, gen::ExternalDecayDriver> HydjetGeneratorFilter;
 }

GeneratorInterface/HydjetInterface/interface/HydjetHadronizer.h

@@ -1,19 +1,24 @@
 #ifndef HydjetHadronizer_h
 #define HydjetHadronizer_h
 
-/** \class HydjetHadronizer
-*
-* Generates HYDJET ==> HepMC events
-*
-* Yetkin Yilmaz
-*   for the Generator Interface. May 2009
-*********************************************/
+/**
+   \class HydjetHadronizer
+   \brief Interface to the HYDJET generator (since core v. 1.9.1), produces HepMC events
+   \version 2.0
+   \authors Camelia Mironov, Yetkin Yilmaz, Andrey Belyaev
+*/
 
 #define PYCOMP pycomp_
 
 #include "FWCore/ParameterSet/interface/ParameterSet.h"
 #include "GeneratorInterface/Core/interface/BaseHadronizer.h"
 
+#include "FWCore/Utilities/interface/EDGetToken.h"
+#include "FWCore/Framework/interface/EDProducer.h"
+#include "FWCore/Framework/interface/ConsumesCollector.h"
+#include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
+#include "SimDataFormats/PileupSummaryInfo/interface/PileupMixingContent.h"
+
 #include <map>
 #include <string>
 #include <vector>
@@ -27,14 +32,15 @@ namespace HepMC {
   class GenEvent;
   class GenParticle;
   class GenVertex;
+  class FourVector;
 }  // namespace HepMC
 
 namespace gen {
   class Pythia6Service;
 
   class HydjetHadronizer : public BaseHadronizer {
   public:
-    HydjetHadronizer(const edm::ParameterSet&);
+    HydjetHadronizer(const edm::ParameterSet&, edm::ConsumesCollector&&);
     ~HydjetHadronizer() override;
 
     bool generatePartonsAndHadronize();
@@ -68,68 +74,69 @@ namespace gen {
 
     HepMC::GenEvent* evt;
     edm::ParameterSet pset_;
-    double abeamtarget_;                  // beam/target atomic mass number
-    int angularspecselector_;             // angular emitted gluon spectrum selection
-    double bfixed_;                       // fixed impact param (fm); valid only if cflag_=0
-    double bmax_;                         // max impact param;
-                                          // units of nucl radius
-    double bmin_;                         // min impact param;
-                                          // units of nucl radius
-    int cflag_;                           // centrality flag
-                                          // =  0 fixed impact param,
-                                          // <> 0 between bmin and bmax
-    bool embedding_;                      // Switch for embedding mode
-    double comenergy;                     // collision energy
-    bool doradiativeenloss_;              //! DEFAULT = true
-    bool docollisionalenloss_;            //! DEFAULT = true
-    double fracsoftmult_;                 // fraction of soft hydro induced hadronic multiplicity
-                                          // proportional to no of nucleon participants
-                                          // (1-fracsoftmult_)--- fraction of soft
-                                          // multiplicity proportional to the numebr
-                                          // of nucleon-nucleon binary collisions
-                                          // DEFAULT=1., allowed range [0.01,1]
-    double hadfreeztemp_;                 // hadron freez-out temperature
-                                          // DEFAULT=0.14MeV, allowed ranges [0.08,0.2]MeV
-    std::string hymode_;                  // Hydjet running mode
-    unsigned int maxEventsToPrint_;       // Events to print if verbosity
-    double maxlongy_;                     // max longitudinal collective rapidity:
-                                          // controls width of eta-spectra
-                                          // DEFAULT=4, allowed range [0.01,7.0]
-    double maxtrany_;                     // max transverse collective rapidity:
-                                          // controls slope of low-pt spectra
-                                          // DEFAULT=1.5, allowed range [0.01,3.0]
-    int nsub_;                            // number of sub-events
-    int nhard_;                           // multiplicity of PYTHIA(+PYQUEN)-induced particles in event
-    int nmultiplicity_;                   // mean soft multiplicity in central PbPb
-                                          // automatically calculated for other centralitie and beams
-    int nsoft_;                           // multiplicity of HYDRO-induced particles in event
-    unsigned int nquarkflavor_;           //! number of active quark flavors in qgp
-                                          //! DEFAULT=0; allowed values: 0,1,2,3.
-    unsigned int pythiaPylistVerbosity_;  // pythia verbosity; def=1
-    double qgpt0_;                        // initial temperature of QGP
-                                          // DEFAULT = 1GeV; allowed range [0.2,2.0]GeV;
-    double qgptau0_;                      // proper time of QGP formation
-                                          // DEFAULT = 0.1 fm/c; allowed range [0.01,10.0]fm/
-
-    double phi0_;  // Event plane angle
+    double abeamtarget_;                  ///< beam/target atomic mass number
+    int angularspecselector_;             ///< angular emitted gluon spectrum selection
+    double bfixed_;                       ///< fixed impact param (fm); valid only if cflag_=0
+    double bmax_;                         ///< max impact param;
+                                          ///< units of nucl radius
+    double bmin_;                         ///< min impact param;
+                                          ///< units of nucl radius
+    int cflag_;                           ///< centrality flag
+                                          ///< =  0 fixed impact param,
+                                          ///< <> 0 between bmin and bmax
+    bool embedding_;                      ///< Switch for embedding mode
+    double comenergy;                     ///< collision energy
+    bool doradiativeenloss_;              ///< DEFAULT = true
+    bool docollisionalenloss_;            ///< DEFAULT = true
+    double fracsoftmult_;                 ///< fraction of soft hydro induced hadronic multiplicity
+                                          ///< proportional to no of nucleon participants
+                                          ///< (1-fracsoftmult_)--- fraction of soft
+                                          ///< multiplicity proportional to the numebr
+                                          ///< of nucleon-nucleon binary collisions
+                                          ///< DEFAULT=1., allowed range [0.01,1]
+    double hadfreeztemp_;                 ///< hadron freez-out temperature
+                                          ///< DEFAULT=0.14MeV, allowed ranges [0.08,0.2]MeV
+    std::string hymode_;                  ///< Hydjet running mode
+    unsigned int maxEventsToPrint_;       ///< Events to print if verbosity
+    double maxlongy_;                     ///< max longitudinal collective rapidity:
+                                          ///< controls width of eta-spectra
+                                          ///< DEFAULT=4, allowed range [0.01,7.0]
+    double maxtrany_;                     ///< max transverse collective rapidity:
+                                          ///< controls slope of low-pt spectra
+                                          ///< DEFAULT=1.5, allowed range [0.01,3.0]
+    int nsub_;                            ///< number of sub-events
+    int nhard_;                           ///< multiplicity of PYTHIA(+PYQUEN)-induced particles in event
+    int nmultiplicity_;                   ///< mean soft multiplicity in central PbPb
+                                          ///< automatically calculated for other centralitie and beams
+    int nsoft_;                           ///< multiplicity of HYDRO-induced particles in event
+    unsigned int nquarkflavor_;           ///< number of active quark flavors in qgp
+                                          ///< DEFAULT=0; allowed values: 0,1,2,3.
+    unsigned int pythiaPylistVerbosity_;  ///< pythia verbosity; def=1
+    double qgpt0_;                        ///< initial temperature of QGP
+                                          ///< DEFAULT = 1GeV; allowed range [0.2,2.0]GeV;
+    double qgptau0_;                      ///< proper time of QGP formation
+                                          ///< DEFAULT = 0.1 fm/c; allowed range [0.01,10.0]fm/
+
+    double phi0_;  ///< Event plane angle
     double sinphi0_;
     double cosphi0_;
-    bool rotate_;  // Switch to rotate event plane
+    bool rotate_;  ///< Switch to rotate event plane
 
-    unsigned int shadowingswitch_;  // shadowing switcher
-                                    // 1-ON, 0-OFF
-    double signn_;                  // inelastic nucleon nucleon cross section [mb]
-                                    // DEFAULT= 58 mb
+    unsigned int shadowingswitch_;   ///< shadowing switcher
+                                     ///< 1-ON, 0-OFF
+    double signn_;                   ///< inelastic nucleon nucleon cross section [mb]
+                                     ///< DEFAULT= 58 mb
+    HepMC::FourVector* fVertex_;     ///< Event signal vertex
+    std::vector<double> signalVtx_;  ///< Pset double vector to set event signal vertex
 
     Pythia6Service* pythia6Service_;
-    edm::InputTag src_;
+    edm::EDGetTokenT<CrossingFrame<edm::HepMCProduct> > src_;
   };
 
-  double HydjetHadronizer::nuclear_radius() const {
-    // Return the nuclear radius derived from the
-    // beam/target atomic mass number.
-
-    return 1.15 * pow((double)abeamtarget_, 1. / 3.);
-  }
+  /**
+    Return the nuclear radius derived from the
+    beam/target atomic mass number.
+ */
+  double HydjetHadronizer::nuclear_radius() const { return 1.15 * pow((double)abeamtarget_, 1. / 3.); }
 }  // namespace gen
 #endif

GeneratorInterface/HydjetInterface/interface/HydjetWrapper.h

@@ -1,16 +1,11 @@
 #ifndef GeneratorInterface_HydjetInterface_HydjetWrapper
 #define GeneratorInterface_HydjetInterface_HydjetWrapper
 
-//
-//
-
-/*
- *
- * Wrapper for FORTRAN version of HYDJET
- *
- * Camelia Mironov
- *
- */
+/**
+   \brief Wrapper for FORTRAN version of HYDJET
+   \version 2.0
+   \authors Camelia Mironov
+*/
 
 extern "C" {
 void hyinit_(double& energy, double& a, int& ifb1, double& bmin, double& bmax, double& bfix1, int& nh1);
@@ -20,10 +15,15 @@ void hyinit_(double& energy, double& a, int& ifb1, double& bmin, double& bmax, d
 #define _MAXMULsize_ 150000
 
 extern "C" {
-void hyevnt_();
+void hyevnt_(double& bfix1);
 }
 #define HYEVNT hyevnt_
 
+extern "C" {
+void hyjver_(int&, int&, int&, int&);
+}
+#define HYJVER hyjver_
+
 extern "C" {
 extern struct { double psi3; } hypsi3_;
 }

GeneratorInterface/HydjetInterface/python/hydjetDefaultParameters_cff.py

@@ -16,7 +16,7 @@
                                    comEnergy = cms.double(5020.0)
                                    )
 
-qgpParameters = cms.PSet(qgpInitialTemperature = cms.double(1.1), ## initial temperature of QGP; allowed range [0.2,2.0]GeV;
+qgpParameters = cms.PSet(qgpInitialTemperature = cms.double(1.), ## initial temperature of QGP; allowed range [0.2,2.0]GeV;
                          qgpProperTimeFormation = cms.double(0.1), ## proper time of QGP formation; allowed range [0.01,10.0]fm/c;
                          hadronFreezoutTemperature = cms.double(0.125),
                          doRadiativeEnLoss = cms.bool(True), ## if true, perform partonic radiative en loss
@@ -27,20 +27,20 @@
 
 hydjetParameters = cms.PSet(sigmaInelNN = cms.double(70),
                             shadowingSwitch = cms.int32(1),
-                            nMultiplicity = cms.int32(25000),
+                            nMultiplicity = cms.int32(18545),
                             fracSoftMultiplicity = cms.double(1.),
-                            maxLongitudinalRapidity = cms.double(4.5),
-                            maxTransverseRapidity = cms.double(1.25),
+                            maxLongitudinalRapidity = cms.double(3.75),
+                            maxTransverseRapidity = cms.double(1.160),
                             rotateEventPlane = cms.bool(True),
                             allowEmptyEvents = cms.bool(False),
-                            angularSpectrumSelector = cms.int32(1), ## angular emitted gluon spectrum :
+                            angularSpectrumSelector = cms.int32(0), ## angular emitted gluon spectrum :
                             embeddingMode = cms.bool(False)
                             )
 
 pyquenPythiaDefaultBlock = cms.PSet(
     pythiaUESettingsBlock,
     hydjetPythiaDefault = cms.vstring('MSEL=0   ! user processes',
-                                      'CKIN(3)=10.',# ! ptMin
+                                      'CKIN(3)=9.2',# ! ptMin
                                       'MSTP(81)=1'
                                       ),
     ppDefault = cms.vstring('MSEL=1   ! QCD hight pT processes (only jets)',