Pythia8Interface/plugins/Py8MassGun.cc

0001
0002 #include <memory>
0003
0004 #include "GeneratorInterface/Core/interface/GeneratorFilter.h"
0005 #include "GeneratorInterface/ExternalDecays/interface/ExternalDecayDriver.h"
0006
0007 #include "GeneratorInterface/Pythia8Interface/interface/Py8GunBase.h"
0008
0009 namespace gen {
0010
0011   class Py8MassGun : public Py8GunBase {
0012   public:
0013     Py8MassGun(edm::ParameterSet const&);
0014     ~Py8MassGun() override {}
0015
0016     bool generatePartonsAndHadronize() override;
0017     const char* classname() const override;
0018
0019   private:
0020     // PtGun particle(s) characteristics
0021     double fMinEta;
0022     double fMaxEta;
0023     double fMinP;
0024     double fMaxP;
0025     double fMinPt;
0026     double fMaxPt;
0027     double fMinM;
0028     double fMaxM;
0029     int fMomMode;
0030   };
0031
0032   // implementation
0033   //
0034   Py8MassGun::Py8MassGun(edm::ParameterSet const& ps) : Py8GunBase(ps) {
0035     // ParameterSet defpset ;
0036     edm::ParameterSet pgun_params = ps.getParameter<edm::ParameterSet>("PGunParameters");  // , defpset ) ;
0037     fMinEta = pgun_params.getParameter<double>("MinEta");                                  // ,-2.2);
0038     fMaxEta = pgun_params.getParameter<double>("MaxEta");                                  // , 2.2);
0039     fMinP = pgun_params.getParameter<double>("MinP");                                      // ,  0.);
0040     fMaxP = pgun_params.getParameter<double>("MaxP");                                      // ,  0.);
0041     fMinPt = pgun_params.getParameter<double>("MinPt");                                    // ,  0.);
0042     fMaxPt = pgun_params.getParameter<double>("MaxPt");                                    // ,  0.);
0043     fMinM = pgun_params.getParameter<double>("MinM");                                      // ,  0.);
0044     fMaxM = pgun_params.getParameter<double>("MaxM");                                      // ,  0.);
0045     fMomMode = pgun_params.getParameter<int>("MomMode");                                   // ,  1);
0046   }
0047
0048   bool Py8MassGun::generatePartonsAndHadronize() {
0049     fMasterGen->event.reset();
0050     size_t pSize = fPartIDs.size();
0051     if (pSize > 2)
0052       return false;
0053
0054     int NTotParticles = fPartIDs.size();
0055
0056     // energy below is dummy, it is not used
0057     (fMasterGen->event).append(990, -11, 0, 0, 2, 1 + NTotParticles, 0, 0, 0., 0., 0., 15000., 15000.);
0058
0059     // Pick a flat mass range
0060     double phi, eta, the, ee, pp;
0061     double m0 = (fMaxM - fMinM) * randomEngine().flat() + fMinM;
0062     // Global eta
0063     eta = (fMaxEta - fMinEta) * randomEngine().flat() + fMinEta;
0064
0065     if (pSize == 2) {
0066       // Masses.
0067       double m1 = fMasterGen->particleData.m0(fPartIDs[0]);
0068       double m2 = fMasterGen->particleData.m0(fPartIDs[1]);
0069
0070       // Energies and absolute momentum in the rest frame.
0071       if (m1 + m2 > m0)
0072         return false;
0073       double e1 = 0.5 * (m0 * m0 + m1 * m1 - m2 * m2) / m0;
0074       double e2 = 0.5 * (m0 * m0 + m2 * m2 - m1 * m1) / m0;
0075       double pAbs = 0.5 * sqrt((m0 - m1 - m2) * (m0 + m1 + m2) * (m0 + m1 - m2) * (m0 - m1 + m2)) / m0;
0076       // Isotropic angles in rest frame give three-momentum.
0077       double cosTheta = 2. * randomEngine().flat() - 1.;
0078       double sinTheta = sqrt(1. - cosTheta * cosTheta);
0079       phi = 2. * M_PI * randomEngine().flat();
0080
0081       double pX = pAbs * sinTheta * cos(phi);
0082       double pY = pAbs * sinTheta * sin(phi);
0083       double pZ = pAbs * cosTheta;
0084
0085       (fMasterGen->event).append(fPartIDs[0], 1, 1, 0, 0, 0, 0, 0, pX, pY, pZ, e1, m1);
0086       (fMasterGen->event).append(fPartIDs[1], 1, 1, 0, 0, 0, 0, 0, -pX, -pY, -pZ, e2, m2);
0087     } else {
0088       (fMasterGen->event).append(fPartIDs[0], 1, 1, 0, 0, 0, 0, 0, 0.0, 0.0, 0.0, m0, m0);
0089     }
0090
0091     //now the boost (from input params)
0092     if (fMomMode == 0) {
0093       pp = (fMaxP - fMinP) * randomEngine().flat() + fMinP;
0094     } else {
0095       double pT = (fMaxPt - fMinPt) * randomEngine().flat() + fMinPt;
0096       pp = pT * cosh(eta);
0097     }
0098     ee = sqrt(m0 * m0 + pp * pp);
0099
0100     //the boost direction (from input params)
0101     //
0102     phi = (fMaxPhi - fMinPhi) * randomEngine().flat() + fMinPhi;
0103     the = 2. * atan(exp(-eta));
0104
0105     double betaX = pp / ee * std::sin(the) * std::cos(phi);
0106     double betaY = pp / ee * std::sin(the) * std::sin(phi);
0107     double betaZ = pp / ee * std::cos(the);
0108
0109     // boost all particles
0110     //
0111     (fMasterGen->event).bst(betaX, betaY, betaZ);
0112
0113     if (!fMasterGen->next())
0114       return false;
0115
0116     event() = std::make_unique<HepMC::GenEvent>();
0117     return toHepMC.fill_next_event(fMasterGen->event, event().get());
0118   }
0119
0120   const char* Py8MassGun::classname() const { return "Py8MassGun"; }
0121
0122   typedef edm::GeneratorFilter<gen::Py8MassGun, gen::ExternalDecayDriver> Pythia8MassGun;
0123
0124 }  // namespace gen
0125
0126 using gen::Pythia8MassGun;
0127 DEFINE_FWK_MODULE(Pythia8MassGun);