Pythia8Interface/plugins/Py8JetGun.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 Py8JetGun : public Py8GunBase {
0012   public:
0013     Py8JetGun(edm::ParameterSet const&);
0014     ~Py8JetGun() 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 fMinE;
0026     double fMaxE;
0027   };
0028
0029   // implementation
0030   //
0031   Py8JetGun::Py8JetGun(edm::ParameterSet const& ps) : Py8GunBase(ps) {
0032     // ParameterSet defpset ;
0033     edm::ParameterSet pgun_params = ps.getParameter<edm::ParameterSet>("PGunParameters");  // , defpset ) ;
0034     fMinEta = pgun_params.getParameter<double>("MinEta");                                  // ,-2.2);
0035     fMaxEta = pgun_params.getParameter<double>("MaxEta");                                  // , 2.2);
0036     fMinP = pgun_params.getParameter<double>("MinP");                                      // ,  0.);
0037     fMaxP = pgun_params.getParameter<double>("MaxP");                                      // ,  0.);
0038     fMinE = pgun_params.getParameter<double>("MinE");                                      // ,  0.);
0039     fMaxE = pgun_params.getParameter<double>("MaxE");                                      // ,  0.);
0040   }
0041
0042   bool Py8JetGun::generatePartonsAndHadronize() {
0043     fMasterGen->event.reset();
0044
0045     int NTotParticles = fPartIDs.size();
0046
0047     // energy below is dummy, it is not used
0048     (fMasterGen->event).append(990, -11, 0, 0, 2, 1 + NTotParticles, 0, 0, 0., 0., 0., 15000., 15000.);
0049
0050     double totPx = 0.;
0051     double totPy = 0.;
0052     double totPz = 0.;
0053     double totE = 0.;
0054     double totM = 0.;
0055     double phi, eta, the, ee, pp;
0056
0057     for (size_t i = 0; i < fPartIDs.size(); i++) {
0058       int particleID = fPartIDs[i];  // this is PDG
0059
0060       phi = 2. * M_PI * randomEngine().flat();
0061       the = acos(-1. + 2. * randomEngine().flat());
0062
0063       // from input
0064       //
0065       ee = (fMaxE - fMinE) * randomEngine().flat() + fMinE;
0066
0067       double mass = (fMasterGen->particleData).m0(particleID);
0068
0069       pp = sqrt(ee * ee - mass * mass);
0070
0071       double px = pp * sin(the) * cos(phi);
0072       double py = pp * sin(the) * sin(phi);
0073       double pz = pp * cos(the);
0074
0075       if (!((fMasterGen->particleData).isParticle(particleID))) {
0076         particleID = std::fabs(particleID);
0077       }
0078       (fMasterGen->event).append(particleID, 1, 1, 0, 0, 0, 0, 0, px, py, pz, ee, mass);
0079       int eventSize = (fMasterGen->event).size() - 1;
0080       // -log(flat) = exponential distribution
0081       double tauTmp = -(fMasterGen->event)[eventSize].tau0() * log(randomEngine().flat());
0082       (fMasterGen->event)[eventSize].tau(tauTmp);
0083
0084       // values for computing total mass
0085       //
0086       totPx += px;
0087       totPy += py;
0088       totPz += pz;
0089       totE += ee;
0090     }
0091
0092     totM = sqrt(totE * totE - (totPx * totPx + totPy * totPy + totPz * totPz));
0093
0094     //now the boost (from input params)
0095     //
0096     pp = (fMaxP - fMinP) * randomEngine().flat() + fMinP;
0097     ee = sqrt(totM * totM + pp * pp);
0098
0099     //the boost direction (from input params)
0100     //
0101     phi = (fMaxPhi - fMinPhi) * randomEngine().flat() + fMinPhi;
0102     eta = (fMaxEta - fMinEta) * randomEngine().flat() + fMinEta;
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* Py8JetGun::classname() const { return "Py8JetGun"; }
0121
0122   typedef edm::GeneratorFilter<gen::Py8JetGun, gen::ExternalDecayDriver> Pythia8JetGun;
0123
0124 }  // namespace gen
0125
0126 using gen::Pythia8JetGun;
0127 DEFINE_FWK_MODULE(Pythia8JetGun);