Pythia8Interface/src/SuepShower.cc

0001 // This is file contains the code to generate a dark sector shower in a strongly coupled, quasi-conformal hidden valley, often
0002 // referred to as "soft unclustered energy patterns (SUEP)" or "softbomb" events. The shower is generated in its rest frame and
0003 // for a realistic simulation this class needs to be interfaced with an event generator such as madgraph, pythia or herwig.
0004 //
0005 // The algorithm relies on arXiv:1305.5226. See arXiv:1612.00850 for a description of the model.
0006 // Please cite both papers when using this code.
0007 //
0008 //
0009 // Applicability:
0010 // The ratio of the parameters m and T (m/T) should be an O(1) number. For m/T>>1 and m/T<<1 the theoretical description of the shower is not valid.
0011 // The mass of the scalar which initiates the shower should be much larger than the mass of the mesons, in other words M>>m,T, by at least an order of magnitude.
0012 //
0013 // Written by Simon Knapen on 12/22/2019, following 1205.5226
0014 // Adapted by Carlos Erice for CMSSW
0015 #include "GeneratorInterface/Pythia8Interface/interface/SuepShower.h"
0016
0017 // constructor
0018 SuepShower::SuepShower(double mass, double temperature, Pythia8::Rndm* rndmPtr) {
0019   edm::LogInfo("Pythia8Interface") << "Creating a SuepShower module";
0020   // Parameter setup
0021   darkmeson_mass_ = mass;
0022   fRndmPtr_ = rndmPtr;
0023   // No need to save temperature, as everything depends on it through m/T
0024   mass_over_T_ = darkmeson_mass_ / temperature;
0025
0026   // 128 bit numerical precision, i.e. take machine precision
0027   tolerance_ = boost::math::tools::eps_tolerance<double>(128);
0028
0029   // Median momentum
0030   p_m_ = sqrt(2 / (mass_over_T_ * mass_over_T_) * (1 + sqrt(1 + mass_over_T_ * mass_over_T_)));
0031   double pmax = sqrt(2 + 2 * sqrt(1 + mass_over_T_ * mass_over_T_)) / mass_over_T_;
0032   // Find the two cases where f(p)/f(pmax) = 1/e, given MB distribution, one appears at each side of pmax always
0033   p_plus_ = (boost::math::tools::bisect(
0034                  boost::bind(&SuepShower::logTestFunction, this, boost::placeholders::_1), pmax, 50.0, tolerance_))
0035                 .first;  // first root
0036   p_minus_ = (boost::math::tools::bisect(
0037                   boost::bind(&SuepShower::logTestFunction, this, boost::placeholders::_1), 0.0, pmax, tolerance_))
0038                  .first;  // second root
0039   // Define the auxiliar quantities for the random generation of momenta
0040   lambda_plus_ = -fMaxwellBoltzmann(p_plus_) / fMaxwellBoltzmannPrime(p_plus_);
0041   lambda_minus_ = fMaxwellBoltzmann(p_minus_) / fMaxwellBoltzmannPrime(p_minus_);
0042   q_plus_ = lambda_plus_ / (p_plus_ - p_minus_);
0043   q_minus_ = lambda_minus_ / (p_plus_ - p_minus_);
0044   q_m_ = 1 - (q_plus_ + q_minus_);
0045 }
0046
0047 SuepShower::~SuepShower() {}
0048
0049 // Generate a shower in the rest frame of the mediator
0050 std::vector<Pythia8::Vec4> SuepShower::generateShower(double energy) {
0051   mediator_energy_ = energy;
0052   std::vector<Pythia8::Vec4> shower;
0053   double shower_energy = 0.0;
0054
0055   // Fill up shower record
0056   while (shower_energy < mediator_energy_ || shower.size() < 2) {
0057     shower.push_back(generateFourVector());
0058     shower_energy += (shower.back()).e();
0059   }
0060
0061   // Reballance momenta to ensure conservation
0062   // Correction starts at 0
0063   Pythia8::Vec4 correction = Pythia8::Vec4(0., 0., 0., 0.);
0064   for (const auto& daughter : shower) {
0065     correction = correction + daughter;
0066   }
0067   correction = correction / shower.size();
0068   // We only want to correct momenta first
0069   correction.e(0);
0070   for (auto& daughter : shower) {
0071     daughter = daughter - correction;
0072   }
0073
0074   // With momentum conserved, balance energy. scale is the multiplicative factor needed such that sum_daughters((scale*p)^2+m^2) = E_parent, i.e. energy is conserved
0075   double scale;
0076   double minscale = 0.0;
0077   double maxscale = 2.0;
0078   while (SuepShower::reballanceFunction(minscale, shower) * SuepShower::reballanceFunction(maxscale, shower) > 0) {
0079     minscale = maxscale;
0080     maxscale *= 2;
0081   }
0082
0083   scale =
0084       (boost::math::tools::bisect(boost::bind(&SuepShower::reballanceFunction, this, boost::placeholders::_1, shower),
0085                                   minscale,
0086                                   maxscale,
0087                                   tolerance_))
0088           .first;
0089
0090   for (auto& daughter : shower) {
0091     daughter.px(daughter.px() * scale);
0092     daughter.py(daughter.py() * scale);
0093     daughter.pz(daughter.pz() * scale);
0094     // Force everything to be on-shell
0095     daughter.e(sqrt(daughter.pAbs2() + darkmeson_mass_ * darkmeson_mass_));
0096   }
0097   return shower;
0098 }
0099
0100 //////// Private Methods ////////
0101 // Maxwell-boltzman distribution, slightly massaged
0102 const double SuepShower::fMaxwellBoltzmann(double p) {
0103   return p * p * exp(-mass_over_T_ * p * p / (1 + sqrt(1 + p * p)));
0104 }
0105
0106 // Derivative of maxwell-boltzmann
0107 const double SuepShower::fMaxwellBoltzmannPrime(double p) {
0108   return exp(-mass_over_T_ * p * p / (1 + sqrt(1 + p * p))) * p * (2 - mass_over_T_ * p * p / sqrt(1 + p * p));
0109 }
0110
0111 // Test function to be solved for p_plus_ and p_minus_
0112 const double SuepShower::logTestFunction(double p) { return log(fMaxwellBoltzmann(p) / fMaxwellBoltzmann(p_m_)) + 1.0; }
0113
0114 // Generate one random 4 std::vector from the thermal distribution
0115 const Pythia8::Vec4 SuepShower::generateFourVector() {
0116   double en, phi, theta, momentum;  //kinematic variables of the 4 std::vector
0117
0118   // First do momentum, following naming of arxiv:1305.5226
0119   double U, V, X, Y, E;
0120   int i = 0;
0121   while (i < 100) {
0122     // Very rarely (<1e-5 events) takes more than one loop to converge, set to 100 for safety
0123     U = fRndmPtr_->flat();
0124     V = fRndmPtr_->flat();
0125
0126     if (U < q_m_) {
0127       Y = U / q_m_;
0128       X = (1 - Y) * (p_minus_ + lambda_minus_) + Y * (p_plus_ - lambda_plus_);
0129       if (V < fMaxwellBoltzmann(X) / fMaxwellBoltzmann(p_m_) && X > 0) {
0130         break;
0131       }
0132     } else {
0133       if (U < q_m_ + q_plus_) {
0134         E = -log((U - q_m_) / q_plus_);
0135         X = p_plus_ - lambda_plus_ * (1 - E);
0136         if (V < exp(E) * fMaxwellBoltzmann(X) / fMaxwellBoltzmann(p_m_) && X > 0) {
0137           break;
0138         }
0139       } else {
0140         E = -log((U - (q_m_ + q_plus_)) / q_minus_);
0141         X = p_minus_ + lambda_minus_ * (1 - E);
0142         if (V < exp(E) * fMaxwellBoltzmann(X) / fMaxwellBoltzmann(p_m_) && X > 0) {
0143           break;
0144         }
0145       }
0146     }
0147   }
0148   // X is the dimensionless momentum, p/m
0149   momentum = X * darkmeson_mass_;
0150
0151   // now do the angles, isotropic
0152   phi = 2.0 * M_PI * (fRndmPtr_->flat());
0153   theta = acos(2.0 * (fRndmPtr_->flat()) - 1.0);
0154
0155   // compose the 4 std::vector
0156   en = sqrt(momentum * momentum + darkmeson_mass_ * darkmeson_mass_);
0157   Pythia8::Vec4 daughterFourMomentum =
0158       Pythia8::Vec4(momentum * cos(phi) * sin(theta), momentum * sin(phi) * sin(theta), momentum * sin(theta), en);
0159   return daughterFourMomentum;
0160 }
0161
0162 // Auxiliary function, to be solved in order to impose energy conservation
0163 // To ballance energy, we solve for "scale" by demanding that this function vanishes, i.e. that shower energy and mediator energy are equal
0164 // By rescaling the momentum rather than the energy, avoid having to do these annoying rotations from the previous version
0165 const double SuepShower::reballanceFunction(double scale, const std::vector<Pythia8::Vec4>& shower) {
0166   double showerEnergy = 0.0;
0167   for (const auto& daughter : shower) {
0168     showerEnergy += sqrt(scale * scale * daughter.pAbs2() + darkmeson_mass_ * darkmeson_mass_);
0169   }
0170   return showerEnergy - mediator_energy_;
0171 }