Pythia8Interface/plugins/Py8PtExpGun.cc

0001 #include <memory>
0002
0003 #include "GeneratorInterface/Core/interface/GeneratorFilter.h"
0004 #include "GeneratorInterface/ExternalDecays/interface/ExternalDecayDriver.h"
0005
0006 #include "GeneratorInterface/Pythia8Interface/interface/Py8GunBase.h"
0007
0008 namespace gen {
0009
0010   class Py8PtExpGun : public Py8GunBase {
0011   public:
0012     Py8PtExpGun(edm::ParameterSet const&);
0013     ~Py8PtExpGun() override {}
0014
0015     bool generatePartonsAndHadronize() override;
0016     const char* classname() const override;
0017
0018   private:
0019     // PtExpGun particle(s) characteristics
0020     double fMinEta;
0021     double fMaxEta;
0022     double fMinPt;
0023     double fMaxPt;
0024     bool fAddAntiParticle;
0025   };
0026
0027   // implementation
0028   //
0029   Py8PtExpGun::Py8PtExpGun(edm::ParameterSet const& ps) : Py8GunBase(ps) {
0030     // ParameterSet defpset ;
0031     edm::ParameterSet pgun_params = ps.getParameter<edm::ParameterSet>("PGunParameters");  // , defpset ) ;
0032     fMinEta = pgun_params.getParameter<double>("MinEta");                                  // ,-2.2);
0033     fMaxEta = pgun_params.getParameter<double>("MaxEta");                                  // , 2.2);
0034     fMinPt = pgun_params.getParameter<double>("MinPt");                                    // ,  0.);
0035     fMaxPt = pgun_params.getParameter<double>("MaxPt");                                    // ,  0.);
0036     fAddAntiParticle = pgun_params.getParameter<bool>("AddAntiParticle");                  //, false) ;
0037   }
0038
0039   bool Py8PtExpGun::generatePartonsAndHadronize() {
0040     fMasterGen->event.reset();
0041
0042     int NTotParticles = fPartIDs.size();
0043     if (fAddAntiParticle)
0044       NTotParticles *= 2;
0045
0046     // energy below is dummy, it is not used
0047     (fMasterGen->event).append(990, -11, 0, 0, 2, 1 + NTotParticles, 0, 0, 0., 0., 0., 15000., 15000.);
0048
0049     for (size_t i = 0; i < fPartIDs.size(); i++) {
0050       int particleID = fPartIDs[i];  // this is PDG - need to convert to Py8 ???
0051
0052       double phi = (fMaxPhi - fMinPhi) * randomEngine().flat() + fMinPhi;
0053       double eta = (fMaxEta - fMinEta) * randomEngine().flat() + fMinEta;
0054       double the = 2. * atan(exp(-eta));
0055
0056       //-log(flat) = exponential distribution
0057       //  need the /10.0 and the min with 1.0 to make sure pt doesn't go too high
0058       //      10.0 chosen to give last pt bin (overflow bin) a reasonable (not unnaturally high) content
0059       double pt = (std::min(-1 / 10.0 * log(randomEngine().flat()), 1.0)) * (fMaxPt - fMinPt) + fMinPt;
0060
0061       double mass = (fMasterGen->particleData).m0(particleID);
0062
0063       double pp = pt / sin(the);  // sqrt( ee*ee - mass*mass );
0064       double ee = sqrt(pp * pp + mass * mass);
0065
0066       double px = pt * cos(phi);
0067       double py = pt * sin(phi);
0068       double pz = pp * cos(the);
0069
0070       if (!((fMasterGen->particleData).isParticle(particleID))) {
0071         particleID = std::abs(particleID);
0072       }
0073       if (1 <= std::abs(particleID) && std::abs(particleID) <= 6)  // quarks
0074         (fMasterGen->event).append(particleID, 23, 1, 0, 0, 0, 101, 0, px, py, pz, ee, mass);
0075       else if (std::abs(particleID) == 21)  // gluons
0076         (fMasterGen->event).append(21, 23, 1, 0, 0, 0, 101, 102, px, py, pz, ee, mass);
0077       // other
0078       else {
0079         (fMasterGen->event).append(particleID, 1, 1, 0, 0, 0, 0, 0, px, py, pz, ee, mass);
0080         int eventSize = (fMasterGen->event).size() - 1;
0081         // -log(flat) = exponential distribution
0082         double tauTmp = -(fMasterGen->event)[eventSize].tau0() * log(randomEngine().flat());
0083         (fMasterGen->event)[eventSize].tau(tauTmp);
0084       }
0085
0086       // Here also need to add anti-particle (if any)
0087       // otherwise just add a 2nd particle of the same type
0088       // (for example, gamma)
0089       //
0090       if (fAddAntiParticle) {
0091         if (1 <= std::abs(particleID) && std::abs(particleID) <= 6) {  // quarks
0092           (fMasterGen->event).append(-particleID, 23, 1, 0, 0, 0, 0, 101, -px, -py, -pz, ee, mass);
0093         } else if (std::abs(particleID) == 21) {  // gluons
0094           (fMasterGen->event).append(21, 23, 1, 0, 0, 0, 102, 101, -px, -py, -pz, ee, mass);
0095         } else {
0096           if ((fMasterGen->particleData).isParticle(-particleID)) {
0097             (fMasterGen->event).append(-particleID, 1, 1, 0, 0, 0, 0, 0, -px, -py, -pz, ee, mass);
0098           } else {
0099             (fMasterGen->event).append(particleID, 1, 1, 0, 0, 0, 0, 0, -px, -py, -pz, ee, mass);
0100           }
0101           int eventSize = (fMasterGen->event).size() - 1;
0102           // -log(flat) = exponential distribution
0103           double tauTmp = -(fMasterGen->event)[eventSize].tau0() * log(randomEngine().flat());
0104           (fMasterGen->event)[eventSize].tau(tauTmp);
0105         }
0106       }
0107     }
0108
0109     if (!fMasterGen->next())
0110       return false;
0111     evtGenDecay();
0112
0113     event() = std::make_unique<HepMC::GenEvent>();
0114     return toHepMC.fill_next_event(fMasterGen->event, event().get());
0115   }
0116
0117   const char* Py8PtExpGun::classname() const { return "Py8PtExpGun"; }
0118
0119   typedef edm::GeneratorFilter<gen::Py8PtExpGun, gen::ExternalDecayDriver> Pythia8PtExpGun;
0120
0121 }  // namespace gen
0122
0123 using gen::Pythia8PtExpGun;
0124 DEFINE_FWK_MODULE(Pythia8PtExpGun);