Pythia8Interface/plugins/Py8PtAndLxyGun.cc

0001 #include <memory>
0002 #include <algorithm>
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 Py8PtAndLxyGun : public Py8GunBase {
0012   public:
0013     Py8PtAndLxyGun(edm::ParameterSet const&);
0014     ~Py8PtAndLxyGun() override {}
0015
0016     bool generatePartonsAndHadronize() override;
0017     const char* classname() const override;
0018
0019   private:
0020     // PtAndLxyGun particle(s) characteristics
0021     double fMinEta;
0022     double fMaxEta;
0023     double fMinPt;
0024     double fMaxPt;
0025     bool fAddAntiParticle;
0026     double fDxyMax;
0027     double fDzMax;
0028     double fLxyMin;
0029     double fLxyMax;
0030     double fLzMax;
0031     double fConeRadius;
0032     double fConeH;
0033     double fDistanceToAPEX;
0034     double fLxyBackFraction;
0035     double fLzOppositeFraction;
0036   };
0037
0038   // implementation
0039   //
0040   Py8PtAndLxyGun::Py8PtAndLxyGun(edm::ParameterSet const& ps) : Py8GunBase(ps) {
0041     edm::ParameterSet pgun_params = ps.getParameter<edm::ParameterSet>("PGunParameters");
0042     fMinEta = pgun_params.getParameter<double>("MinEta");
0043     fMaxEta = pgun_params.getParameter<double>("MaxEta");
0044     fMinPt = pgun_params.getParameter<double>("MinPt");
0045     fMaxPt = pgun_params.getParameter<double>("MaxPt");
0046     fAddAntiParticle = pgun_params.getParameter<bool>("AddAntiParticle");
0047     fDxyMax = pgun_params.getParameter<double>("dxyMax");
0048     fDzMax = pgun_params.getParameter<double>("dzMax");
0049     fLxyMin = pgun_params.getParameter<double>("LxyMin");
0050     fLxyMax = pgun_params.getParameter<double>("LxyMax");
0051     fLzMax = pgun_params.getParameter<double>("LzMax");
0052     fConeRadius = pgun_params.getParameter<double>("ConeRadius");
0053     fConeH = pgun_params.getParameter<double>("ConeH");
0054     fDistanceToAPEX = pgun_params.getParameter<double>("DistanceToAPEX");
0055     fLxyBackFraction = std::clamp(pgun_params.getParameter<double>("LxyBackFraction"), 0., 1.);
0056     fLzOppositeFraction = std::clamp(pgun_params.getParameter<double>("LzOppositeFraction"), 0., 1.);
0057   }
0058
0059   bool Py8PtAndLxyGun::generatePartonsAndHadronize() {
0060     fMasterGen->event.reset();
0061
0062     for (size_t i = 0; i < fPartIDs.size(); i++) {
0063       int particleID = fPartIDs[i];  // this is PDG - need to convert to Py8 ???
0064
0065       double phi = 0;
0066       double dxy = 0;
0067       double pt = 0;
0068       double eta = 0;
0069       double px = 0;
0070       double py = 0;
0071       double pz = 0;
0072       double mass = 0;
0073       double ee = 0;
0074       double vx = 0;
0075       double vy = 0;
0076       double vz = 0;
0077       double lxy = 0;
0078
0079       bool passLoop = false;
0080       while (!passLoop) {
0081         bool passDxy = false;
0082         bool passLz = false;
0083         bool passDz = false;
0084
0085         phi = (fMaxPhi - fMinPhi) * randomEngine().flat() + fMinPhi;
0086         pt = (fMaxPt - fMinPt) * randomEngine().flat() + fMinPt;
0087         px = pt * cos(phi);
0088         py = pt * sin(phi);
0089
0090         lxy = (fLxyMax - fLxyMin) * randomEngine().flat() + fLxyMin;
0091
0092         int sign = 1;
0093         for (int i = 0; i < 10000; i++) {
0094           double vphi = 2 * M_PI * randomEngine().flat();
0095           vx = lxy * cos(vphi);
0096           vy = lxy * sin(vphi);
0097
0098           dxy = -vx * sin(phi) + vy * cos(phi);
0099
0100           sign = 1;
0101           if (fLxyBackFraction > 0 && randomEngine().flat() <= fLxyBackFraction) {
0102             sign = -1;
0103           }
0104           if ((std::abs(dxy) < fDxyMax || fDxyMax < 0) && sign * (vx * px + vy * py) > 0) {
0105             passDxy = true;
0106             break;
0107           }
0108         }
0109
0110         eta = (fMaxEta - fMinEta) * randomEngine().flat() + fMinEta;
0111         double theta = 2. * atan(exp(-eta));
0112
0113         mass = (fMasterGen->particleData).m0(particleID);
0114
0115         double pp = pt / sin(theta);  // sqrt( ee*ee - mass*mass );
0116         ee = sqrt(pp * pp + mass * mass);
0117
0118         pz = pp * cos(theta);
0119
0120         float coneTheta = fConeRadius / fConeH;
0121         for (int j = 0; j < 100; j++) {
0122           vz = fLzMax * randomEngine().flat();  // this is abs(vz)
0123           float v0 = vz - fDistanceToAPEX;
0124           if (v0 <= 0 || lxy * lxy / (coneTheta * coneTheta) > v0 * v0) {
0125             passLz = true;
0126             break;
0127           }
0128         }
0129
0130         if (fLzOppositeFraction > 0 && randomEngine().flat() <= fLzOppositeFraction)
0131           sign *= -1;
0132         if (sign * pz < 0)
0133           vz = -vz;
0134
0135         double dz = vz - (vx * cos(phi) + vy * sin(phi)) / tan(theta);
0136         if (std::abs(dz) < fDzMax || fDzMax < 0) {
0137           passDz = true;
0138         }
0139
0140         passLoop = (passDxy && passLz && passDz);
0141         if (passLoop)
0142           break;
0143       }
0144
0145       float time = sqrt(vx * vx + vy * vy + vz * vz);
0146
0147       if (!((fMasterGen->particleData).isParticle(particleID))) {
0148         particleID = std::abs(particleID);
0149       }
0150       if (1 <= std::abs(particleID) && std::abs(particleID) <= 6)  // quarks
0151         (fMasterGen->event).append(particleID, 23, 101, 0, px, py, pz, ee, mass);
0152       else if (std::abs(particleID) == 21)  // gluons
0153         (fMasterGen->event).append(21, 23, 101, 102, px, py, pz, ee, mass);
0154       // other
0155       else {
0156         (fMasterGen->event).append(particleID, 1, 0, 0, px, py, pz, ee, mass);
0157         int eventSize = (fMasterGen->event).size() - 1;
0158         // -log(flat) = exponential distribution
0159         double tauTmp = -(fMasterGen->event)[eventSize].tau0() * log(randomEngine().flat());
0160         (fMasterGen->event)[eventSize].tau(tauTmp);
0161       }
0162       (fMasterGen->event).back().vProd(vx, vy, vz, time);
0163
0164       // Here also need to add anti-particle (if any)
0165       // otherwise just add a 2nd particle of the same type
0166       // (for example, gamma).
0167       // Added anti-particle has momentum opposite to corresponding
0168       // particle, (px,py,pz)=>(-px,-py,-pz), and production vertex
0169       // symmetric wrt (0,0,0), (vx, vy, vz)=>(-vx, -vy, -vz).
0170       //
0171       if (fAddAntiParticle) {
0172         if (1 <= std::abs(particleID) && std::abs(particleID) <= 6) {  // quarks
0173           (fMasterGen->event).append(-particleID, 23, 0, 101, -px, -py, -pz, ee, mass);
0174         } else if (std::abs(particleID) == 21) {  // gluons
0175           (fMasterGen->event).append(21, 23, 102, 101, -px, -py, -pz, ee, mass);
0176         } else {
0177           if ((fMasterGen->particleData).isParticle(-particleID)) {
0178             (fMasterGen->event).append(-particleID, 1, 0, 0, -px, -py, -pz, ee, mass);
0179           } else {
0180             (fMasterGen->event).append(particleID, 1, 0, 0, -px, -py, -pz, ee, mass);
0181           }
0182           int eventSize = (fMasterGen->event).size() - 1;
0183           // -log(flat) = exponential distribution
0184           double tauTmp = -(fMasterGen->event)[eventSize].tau0() * log(randomEngine().flat());
0185           (fMasterGen->event)[eventSize].tau(tauTmp);
0186         }
0187         (fMasterGen->event).back().vProd(-vx, -vy, -vz, time);
0188       }
0189     }
0190
0191     if (!fMasterGen->next())
0192       return false;
0193     evtGenDecay();
0194
0195     event() = std::make_unique<HepMC::GenEvent>();
0196     return toHepMC.fill_next_event(fMasterGen->event, event().get());
0197   }
0198
0199   const char* Py8PtAndLxyGun::classname() const { return "Py8PtAndLxyGun"; }
0200
0201   typedef edm::GeneratorFilter<gen::Py8PtAndLxyGun, gen::ExternalDecayDriver> Pythia8PtAndLxyGun;
0202
0203 }  // namespace gen
0204
0205 using gen::Pythia8PtAndLxyGun;
0206 DEFINE_FWK_MODULE(Pythia8PtAndLxyGun);