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