Pythia8Interface/plugins/Py8EGun.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 Py8EGun : public Py8GunBase {
0012   public:
0013     Py8EGun(edm::ParameterSet const&);
0014     ~Py8EGun() override {}
0015
0016     bool generatePartonsAndHadronize() override;
0017     const char* classname() const override;
0018
0019   private:
0020     // EGun particle(s) characteristics
0021     double fMinEta;
0022     double fMaxEta;
0023     double fMinE;
0024     double fMaxE;
0025     bool fAddAntiParticle;
0026   };
0027
0028   // implementation
0029   //
0030   Py8EGun::Py8EGun(edm::ParameterSet const& ps) : Py8GunBase(ps) {
0031     // ParameterSet defpset ;
0032     edm::ParameterSet pgun_params = ps.getParameter<edm::ParameterSet>("PGunParameters");  // , defpset ) ;
0033     fMinEta = pgun_params.getParameter<double>("MinEta");                                  // ,-2.2);
0034     fMaxEta = pgun_params.getParameter<double>("MaxEta");                                  // , 2.2);
0035     fMinE = pgun_params.getParameter<double>("MinE");                                      // ,  0.);
0036     fMaxE = pgun_params.getParameter<double>("MaxE");                                      // ,  0.);
0037     fAddAntiParticle = pgun_params.getParameter<bool>("AddAntiParticle");                  //, false) ;
0038   }
0039
0040   bool Py8EGun::generatePartonsAndHadronize() {
0041     fMasterGen->event.reset();
0042
0043     int NTotParticles = fPartIDs.size();
0044     if (fAddAntiParticle)
0045       NTotParticles *= 2;
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     int colorindex = 101;
0051
0052     for (size_t i = 0; i < fPartIDs.size(); i++) {
0053       int particleID = fPartIDs[i];  // this is PDG - need to convert to Py8 ???
0054       if ((std::abs(particleID) <= 6 || particleID == 21) && !(fAddAntiParticle)) {
0055         throw cms::Exception("PythiaError") << "Attempting to generate quarks or gluons without setting "
0056                                                "AddAntiParticle to true. This will not handle color properly."
0057                                             << std::endl;
0058       }
0059
0060       double phi = (fMaxPhi - fMinPhi) * randomEngine().flat() + fMinPhi;
0061       double ee = (fMaxE - fMinE) * randomEngine().flat() + fMinE;
0062       double eta = (fMaxEta - fMinEta) * randomEngine().flat() + fMinEta;
0063       double the = 2. * atan(exp(-eta));
0064
0065       double mass = (fMasterGen->particleData).m0(particleID);
0066
0067       double pp = sqrt(ee * ee - mass * mass);
0068       double px = pp * sin(the) * cos(phi);
0069       double py = pp * sin(the) * sin(phi);
0070       double pz = pp * cos(the);
0071
0072       if (!((fMasterGen->particleData).isParticle(particleID))) {
0073         particleID = std::fabs(particleID);
0074       }
0075       if (1 <= std::abs(particleID) && std::abs(particleID) <= 6) {  // quarks
0076         (fMasterGen->event).append(particleID, 23, 1, 0, 0, 0, colorindex, 0, px, py, pz, ee, mass);
0077         if (!fAddAntiParticle)
0078           colorindex += 1;
0079       } else if (std::abs(particleID) == 21) {  // gluons
0080         (fMasterGen->event).append(21, 23, 1, 0, 0, 0, colorindex, colorindex + 1, px, py, pz, ee, mass);
0081         if (!fAddAntiParticle) {
0082           colorindex += 2;
0083         }
0084       }
0085       // other
0086       else {
0087         (fMasterGen->event).append(particleID, 1, 1, 0, 0, 0, 0, 0, px, py, pz, ee, mass);
0088         int eventSize = (fMasterGen->event).size() - 1;
0089         // -log(flat) = exponential distribution
0090         double tauTmp = -(fMasterGen->event)[eventSize].tau0() * log(randomEngine().flat());
0091         (fMasterGen->event)[eventSize].tau(tauTmp);
0092       }
0093
0094       // Here also need to add anti-particle (if any)
0095       // otherwise just add a 2nd particle of the same type
0096       // (for example, gamma)
0097       //
0098       if (fAddAntiParticle) {
0099         if (1 <= std::abs(particleID) && std::abs(particleID) <= 6) {  // quarks
0100           (fMasterGen->event).append(-particleID, 23, 1, 0, 0, 0, 0, colorindex, -px, -py, -pz, ee, mass);
0101           colorindex += 1;
0102         } else if (std::abs(particleID) == 21) {  // gluons
0103           (fMasterGen->event).append(21, 23, 1, 0, 0, 0, colorindex + 1, colorindex, -px, -py, -pz, ee, mass);
0104           colorindex += 2;
0105         } else {
0106           if ((fMasterGen->particleData).isParticle(-particleID)) {
0107             (fMasterGen->event).append(-particleID, 1, 1, 0, 0, 0, 0, 0, -px, -py, -pz, ee, mass);
0108           } else {
0109             (fMasterGen->event).append(particleID, 1, 1, 0, 0, 0, 0, 0, -px, -py, -pz, ee, mass);
0110           }
0111           int eventSize = (fMasterGen->event).size() - 1;
0112           // -log(flat) = exponential distribution
0113           double tauTmp = -(fMasterGen->event)[eventSize].tau0() * log(randomEngine().flat());
0114           (fMasterGen->event)[eventSize].tau(tauTmp);
0115         }
0116       }
0117     }
0118
0119     if (!fMasterGen->next())
0120       return false;
0121     evtGenDecay();
0122
0123     event() = std::make_unique<HepMC::GenEvent>();
0124     return toHepMC.fill_next_event(fMasterGen->event, event().get());
0125   }
0126
0127   const char* Py8EGun::classname() const { return "Py8EGun"; }
0128
0129   typedef edm::GeneratorFilter<gen::Py8EGun, gen::ExternalDecayDriver> Pythia8EGun;
0130
0131 }  // namespace gen
0132
0133 using gen::Pythia8EGun;
0134 DEFINE_FWK_MODULE(Pythia8EGun);