Flat mass gun

GeneratorInterface/Pythia8Interface/plugins/Py8MassGun.cc

@@ -0,0 +1,120 @@
+
+#include "GeneratorInterface/Core/interface/GeneratorFilter.h"
+#include "GeneratorInterface/ExternalDecays/interface/ExternalDecayDriver.h"
+
+#include "GeneratorInterface/Pythia8Interface/interface/Py8GunBase.h"
+
+namespace gen {
+
+  class Py8MassGun : public Py8GunBase {
+  public:
+    Py8MassGun(edm::ParameterSet const&);
+    ~Py8MassGun() override {}
+
+    bool generatePartonsAndHadronize() override;
+    const char* classname() const override;
+
+  private:
+    // PtGun particle(s) characteristics
+    double fMinEta;
+    double fMaxEta;
+    double fMinP;
+    double fMaxP;
+    double fMinPt;
+    double fMaxPt;
+    double fMinM;
+    double fMaxM;
+    int fMomMode;
+  };
+
+  // implementation
+  //
+  Py8MassGun::Py8MassGun(edm::ParameterSet const& ps) : Py8GunBase(ps) {
+    // ParameterSet defpset ;
+    edm::ParameterSet pgun_params = ps.getParameter<edm::ParameterSet>("PGunParameters");  // , defpset ) ;
+    fMinEta = pgun_params.getParameter<double>("MinEta");                                  // ,-2.2);
+    fMaxEta = pgun_params.getParameter<double>("MaxEta");                                  // , 2.2);
+    fMinP = pgun_params.getParameter<double>("MinP");                                      // ,  0.);
+    fMaxP = pgun_params.getParameter<double>("MaxP");                                      // ,  0.);
+    fMinPt = pgun_params.getParameter<double>("MinPt");                                    // ,  0.);
+    fMaxPt = pgun_params.getParameter<double>("MaxPt");                                    // ,  0.);
+    fMinM = pgun_params.getParameter<double>("MinM");                                      // ,  0.);
+    fMaxM = pgun_params.getParameter<double>("MaxM");                                      // ,  0.);
+    fMomMode = pgun_params.getParameter<int>("MomMode");                                   // ,  1);
+  }
+
+  bool Py8MassGun::generatePartonsAndHadronize() {
+    fMasterGen->event.reset();
+    size_t pSize = fPartIDs.size();
+    if (pSize > 2)
+      return false;
+
+    // Pick a flat mass range
+    double phi, eta, the, ee, pp;
+    double m0 = (fMaxM - fMinM) * randomEngine().flat() + fMinM;
+    // Global eta
+    eta = (fMaxEta - fMinEta) * randomEngine().flat() + fMinEta;
+
+    if (pSize == 2) {
+      // Masses.
+      double m1 = fMasterGen->particleData.m0(fPartIDs[0]);
+      double m2 = fMasterGen->particleData.m0(fPartIDs[1]);
+
+      // Energies and absolute momentum in the rest frame.
+      if (m1 + m2 > m0)
+        return false;
+      double e1 = 0.5 * (m0 * m0 + m1 * m1 - m2 * m2) / m0;
+      double e2 = 0.5 * (m0 * m0 + m2 * m2 - m1 * m1) / m0;
+      double pAbs = 0.5 * sqrt((m0 - m1 - m2) * (m0 + m1 + m2) * (m0 + m1 - m2) * (m0 - m1 + m2)) / m0;
+      // Isotropic angles in rest frame give three-momentum.
+      double cosTheta = 2. * randomEngine().flat() - 1.;
+      double sinTheta = sqrt(1. - cosTheta * cosTheta);
+      phi = 2. * M_PI * randomEngine().flat();
+
+      double pX = pAbs * sinTheta * cos(phi);
+      double pY = pAbs * sinTheta * sin(phi);
+      double pZ = pAbs * cosTheta;
+
+      (fMasterGen->event).append(fPartIDs[0], 1, 0, 0, pX, pY, pZ, e1, m1);
+      (fMasterGen->event).append(fPartIDs[1], 1, 0, 0, -pX, -pY, -pZ, e2, m2);
+    } else {
+      (fMasterGen->event).append(fPartIDs[0], 1, 0, 0, 0.0, 0.0, 0.0, m0, m0);
+    }
+
+    //now the boost (from input params)
+    if (fMomMode == 0) {
+      pp = (fMaxP - fMinP) * randomEngine().flat() + fMinP;
+    } else {
+      double pT = (fMaxPt - fMinPt) * randomEngine().flat() + fMinPt;
+      pp = pT * cosh(eta);
+    }
+    ee = sqrt(m0 * m0 + pp * pp);
+
+    //the boost direction (from input params)
+    //
+    phi = (fMaxPhi - fMinPhi) * randomEngine().flat() + fMinPhi;
+    the = 2. * atan(exp(-eta));
+
+    double betaX = pp / ee * std::sin(the) * std::cos(phi);
+    double betaY = pp / ee * std::sin(the) * std::sin(phi);
+    double betaZ = pp / ee * std::cos(the);
+
+    // boost all particles
+    //
+    (fMasterGen->event).bst(betaX, betaY, betaZ);
+
+    if (!fMasterGen->next())
+      return false;
+
+    event().reset(new HepMC::GenEvent);
+    return toHepMC.fill_next_event(fMasterGen->event, event().get());
+  }
+
+  const char* Py8MassGun::classname() const { return "Py8MassGun"; }
+
+  typedef edm::GeneratorFilter<gen::Py8MassGun, gen::ExternalDecayDriver> Pythia8MassGun;
+
+}  // namespace gen
+
+using gen::Pythia8MassGun;
+DEFINE_FWK_MODULE(Pythia8MassGun);

GeneratorInterface/Pythia8Interface/test/Py8MassGun_cfg.py

@@ -0,0 +1,85 @@
+import FWCore.ParameterSet.Config as cms
+
+from Configuration.Generator.Pythia8CommonSettings_cfi import *
+
+process = cms.Process("TEST")
+
+process.load("FWCore.Framework.test.cmsExceptionsFatal_cff")
+process.load("SimGeneral.HepPDTESSource.pythiapdt_cfi")
+
+process.source = cms.Source("EmptySource")
+
+process.generator = cms.EDFilter("Pythia8MassGun", 
+    maxEventsToPrint = cms.untracked.int32(100),
+    pythiaPylistVerbosity = cms.untracked.int32(1),
+    pythiaHepMCVerbosity = cms.untracked.bool(True),
+
+    PGunParameters = cms.PSet(
+       ParticleID = cms.vint32(999999),
+        # this defines "absolute" energy spead of particles in the jet
+	MinM   = cms.double(8.0),
+	MaxM   = cms.double(15.0),
+	# the following params define the boost
+        MinP   = cms.double(20.0),
+        MaxP   = cms.double(20.0),
+        MomMode = cms.int32(1),
+        MinPt   = cms.double(20.0),
+        MaxPt   = cms.double(20.0),
+        MinPhi = cms.double(-3.14159265359),
+        MaxPhi = cms.double(3.14159265359),
+	MinEta = cms.double(-2.4),
+        MaxEta = cms.double(2.4)
+    ),
+    PythiaParameters = cms.PSet(
+
+    pythia8CommonSettingsBlock,
+
+    processParameters = cms.vstring(
+      '999999:all = GeneralResonance void 1 0 0 500. 1. 0. 0. 0.',
+      '999999:oneChannel = 1 1.00 101 15 -15 15 -15',
+      'Main:timesAllowErrors    = 10000',
+      '15:onMode = off',
+      '15:onIfAll = 211 211 211',
+      '15:onIfAll = 211 211 321',
+      '15:onIfAll = 211 321 321',
+      '15:onIfAll = 321 321 321',
+      '15:onIfAll = 321 321 321',
+      ),
+
+    parameterSets = cms.vstring(
+        'processParameters',
+	'pythia8CommonSettings')
+    )
+)
+
+process.load("FWCore.MessageLogger.MessageLogger_cfi")
+process.MessageLogger = cms.Service("MessageLogger",
+    cout = cms.untracked.PSet(
+        default = cms.untracked.PSet(
+            limit = cms.untracked.int32(2)
+        )
+    ),
+    destinations = cms.untracked.vstring('cout')
+)
+
+process.RandomNumberGeneratorService = cms.Service("RandomNumberGeneratorService",
+    generator = cms.PSet(
+        initialSeed = cms.untracked.uint32(123456789),
+        engineName = cms.untracked.string('HepJamesRandom')
+    )
+)
+
+process.maxEvents = cms.untracked.PSet(
+    input = cms.untracked.int32(100)
+)
+
+process.GEN = cms.OutputModule("PoolOutputModule",
+    fileName = cms.untracked.string('Py8JetGun.root')
+)
+
+
+process.p = cms.Path(process.generator)
+process.outpath = cms.EndPath(process.GEN)
+
+process.schedule = cms.Schedule(process.p, process.outpath)
+